From 1e94a9dfe39187b3eaca31b91a9f83f7b0100aba Mon Sep 17 00:00:00 2001 From: Jan Wilhelm Date: Tue, 7 Dec 2021 13:49:23 +0100 Subject: [PATCH] regularized RI for periodic GW; diagonalization and inversion instead of Cholesky decomposition --- src/fm/cp_cfm_basic_linalg.F | 4 +- src/fm/cp_cfm_diag.F | 2 +- src/input_constants.F | 12 +- src/input_cp2k_dft.F | 31 +- src/input_cp2k_kpoints.F | 41 +++ src/input_cp2k_mp2.F | 126 +++++-- src/kpoint_coulomb_2c.F | 18 +- src/kpoint_methods.F | 39 ++- src/mp2_integrals.F | 43 ++- src/mp2_ri_2c.F | 120 +++---- src/mp2_setup.F | 65 +++- src/mp2_types.F | 14 +- src/qs_band_structure.F | 48 +-- src/qs_gamma2kp.F | 57 ++- src/qs_ks_methods.F | 1 + src/qs_wannier90.F | 28 +- src/rpa_gw.F | 27 +- src/rpa_gw_im_time_util.F | 177 +++++++++- src/rpa_gw_kpoints.F | 330 +++++++----------- src/rpa_im_time.F | 280 +++++++++++++-- src/rpa_kpoints.F | 225 ++++++++---- src/rpa_main.F | 213 +++++++++-- src/rpa_util.F | 54 +-- .../G0W0_H2O_PBE_periodic.inp | 6 +- ...ints_from_Gamma_RI_regularization_1E-3.inp | 83 +++++ tests/QS/regtest-gw-cubic/TEST_FILES | 3 +- tests/QS/regtest-gw/G0W0_H2O_PBE_periodic.inp | 7 +- 27 files changed, 1437 insertions(+), 617 deletions(-) create mode 100644 tests/QS/regtest-gw-cubic/G0W0_kpoints_from_Gamma_RI_regularization_1E-3.inp diff --git a/src/fm/cp_cfm_basic_linalg.F b/src/fm/cp_cfm_basic_linalg.F index ff9ffe29b2..2b3b5df3f1 100644 --- a/src/fm/cp_cfm_basic_linalg.F +++ b/src/fm/cp_cfm_basic_linalg.F @@ -844,6 +844,7 @@ END SUBROUTINE cp_cfm_lu_invert !> \par History !> 05.2002 created [JVdV] !> 12.2002 updated, added n optional parm [fawzi] +!> 09.2021 removed CPASSERT(info == 0) since there is already check of info [Jan Wilhelm] !> \author Joost ! ************************************************************************************************** SUBROUTINE cp_cfm_cholesky_decompose(matrix, n, info_out) @@ -888,9 +889,8 @@ SUBROUTINE cp_cfm_cholesky_decompose(matrix, n, info_out) "Cholesky decompose failed: matrix is not positive definite or ill-conditioned") END IF - CPASSERT(info == 0) - CALL timestop(handle) + END SUBROUTINE cp_cfm_cholesky_decompose ! ************************************************************************************************** diff --git a/src/fm/cp_cfm_diag.F b/src/fm/cp_cfm_diag.F index 1e8c891489..7526e91c84 100644 --- a/src/fm/cp_cfm_diag.F +++ b/src/fm/cp_cfm_diag.F @@ -212,7 +212,7 @@ SUBROUTINE cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, CALL timeset(routineN, handle) - ! Test sizees + ! Test sizes CALL cp_cfm_get_info(amatrix, nrow_global=nao) nmo = SIZE(eigenvalues) ALLOCATE (evals(nao), cevals(nao)) diff --git a/src/input_constants.F b/src/input_constants.F index 72f0a484f9..755a8e2e95 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -1006,17 +1006,19 @@ MODULE input_constants wfc_mm_style_gemm = 11, & wfc_mm_style_syrk = 12 - ! G0W0 parameter + ! GW parameter INTEGER, PARAMETER, PUBLIC :: ri_rpa_g0w0_crossing_z_shot = 1, & ri_rpa_g0w0_crossing_newton = 2, & ri_rpa_g0w0_crossing_bisection = 3, & gw_no_print_exx = 5, & gw_print_exx = 6, & gw_read_exx = 7, & - gw_skip_for_regtest = 8 - - INTEGER, PARAMETER, PUBLIC :: gw_pade_approx = 0, & - gw_two_pole_model = 1 + gw_skip_for_regtest = 8, & + gw_pade_approx = 9, & + gw_two_pole_model = 10, & + kp_weights_W_auto = 11, & + kp_weights_W_uniform = 12, & + kp_weights_W_tailored = 13 ! periodic RESP parameters INTEGER, PARAMETER, PUBLIC :: do_resp_x_dir = 0, & diff --git a/src/input_cp2k_dft.F b/src/input_cp2k_dft.F index 37efd821b3..8b2ee1b00e 100644 --- a/src/input_cp2k_dft.F +++ b/src/input_cp2k_dft.F @@ -110,7 +110,8 @@ MODULE input_cp2k_dft create_ext_vxc_section USE input_cp2k_field, ONLY: create_efield_section,& create_per_efield_section - USE input_cp2k_kpoints, ONLY: create_kpoints_section + USE input_cp2k_kpoints, ONLY: create_kpoint_set_section,& + create_kpoints_section USE input_cp2k_loc, ONLY: create_localize_section,& print_wanniers USE input_cp2k_ls, ONLY: create_ls_scf_section @@ -2340,33 +2341,7 @@ SUBROUTINE create_bandstructure_section(section) CALL keyword_release(keyword) NULLIFY (subsection) - CALL section_create(subsection, __LOCATION__, name="KPOINT_SET", & - description="Specifies a k-point line to be calculated.", & - n_keywords=0, n_subsections=0, repeats=.TRUE.) - ! keywords - NULLIFY (keyword) - CALL keyword_create(keyword, __LOCATION__, name="SPECIAL_POINT", & - description="Name and coordinates of a special k-point", & - usage="SPECIAL_POINT GAMMA 0.0 0.0 0.0", n_var=-1, type_of_var=char_t, repeats=.TRUE.) - CALL section_add_keyword(subsection, keyword) - CALL keyword_release(keyword) - ! - CALL keyword_create(keyword, __LOCATION__, name="NPOINTS", & - description="Number of k-points along the line.", & - default_i_val=0) - CALL section_add_keyword(subsection, keyword) - CALL keyword_release(keyword) - ! - CALL keyword_create(keyword, __LOCATION__, name="UNITS", & - description="Special k-points are defined either in units"// & - " of reciprocal lattice vectors or in Cartesian coordinates in uints of 2Pi/len."// & - " B_VECTOR: in multiples of the reciprocal lattice vectors (b)."// & - " CART_ANGSTROM: In units of 2*Pi/Angstrom."// & - " CART_BOHR: In units of 2*Pi/Bohr.", & - usage="UNITS ", type_of_var=char_t, default_c_val="B_VECTOR") - CALL section_add_keyword(subsection, keyword) - CALL keyword_release(keyword) - ! + CALL create_kpoint_set_section(subsection) CALL section_add_subsection(section, subsection) CALL section_release(subsection) diff --git a/src/input_cp2k_kpoints.F b/src/input_cp2k_kpoints.F index 00e75c9630..c3b7d6b91c 100644 --- a/src/input_cp2k_kpoints.F +++ b/src/input_cp2k_kpoints.F @@ -35,6 +35,7 @@ MODULE input_cp2k_kpoints use_complex_wfn = 100 PUBLIC :: create_kpoints_section + PUBLIC :: create_kpoint_set_section PUBLIC :: use_real_wfn, use_complex_wfn CONTAINS @@ -149,4 +150,44 @@ SUBROUTINE create_kpoints_section(section) END SUBROUTINE create_kpoints_section +! ************************************************************************************************** +!> \brief ... +!> \param section ... +!> \author JGH (moved to a new routine by Jan Wilhelm, 11.2021) +! ************************************************************************************************** + SUBROUTINE create_kpoint_set_section(section) + TYPE(section_type), POINTER :: section + + TYPE(keyword_type), POINTER :: keyword + + CPASSERT(.NOT. ASSOCIATED(section)) + CALL section_create(section, __LOCATION__, name="KPOINT_SET", & + description="Specifies a k-point line to be calculated.", & + n_keywords=0, n_subsections=0, repeats=.TRUE.) + ! keywords + NULLIFY (keyword) + CALL keyword_create(keyword, __LOCATION__, name="SPECIAL_POINT", & + description="Name and coordinates of a special k-point", & + usage="SPECIAL_POINT GAMMA 0.0 0.0 0.0", n_var=-1, type_of_var=char_t, repeats=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + ! + CALL keyword_create(keyword, __LOCATION__, name="NPOINTS", & + description="Number of k-points along the line.", & + default_i_val=0) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + ! + CALL keyword_create(keyword, __LOCATION__, name="UNITS", & + description="Special k-points are defined either in units"// & + " of reciprocal lattice vectors or in Cartesian coordinates in uints of 2Pi/len."// & + " B_VECTOR: in multiples of the reciprocal lattice vectors (b)."// & + " CART_ANGSTROM: In units of 2*Pi/Angstrom."// & + " CART_BOHR: In units of 2*Pi/Bohr.", & + usage="UNITS ", type_of_var=char_t, default_c_val="B_VECTOR") + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + END SUBROUTINE create_kpoint_set_section + END MODULE input_cp2k_kpoints diff --git a/src/input_cp2k_mp2.F b/src/input_cp2k_mp2.F index a63f134afd..9f96e64627 100644 --- a/src/input_cp2k_mp2.F +++ b/src/input_cp2k_mp2.F @@ -28,11 +28,12 @@ MODULE input_cp2k_mp2 do_eri_gpw, do_eri_mme, do_eri_os, do_potential_coulomb, do_potential_id, & do_potential_long, do_potential_short, do_potential_truncated, do_potential_tshpsc, & eri_default, gaussian, gw_no_print_exx, gw_pade_approx, gw_print_exx, gw_read_exx, & - gw_skip_for_regtest, gw_two_pole_model, mp2_method_direct, mp2_method_gpw, & - mp2_method_none, numerical, ri_default, ri_rpa_g0w0_crossing_bisection, & - ri_rpa_g0w0_crossing_newton, ri_rpa_g0w0_crossing_z_shot, wfc_mm_style_gemm, & - wfc_mm_style_syrk + gw_skip_for_regtest, gw_two_pole_model, kp_weights_W_auto, kp_weights_W_tailored, & + kp_weights_W_uniform, mp2_method_direct, mp2_method_gpw, mp2_method_none, numerical, & + ri_default, ri_rpa_g0w0_crossing_bisection, ri_rpa_g0w0_crossing_newton, & + ri_rpa_g0w0_crossing_z_shot, wfc_mm_style_gemm, wfc_mm_style_syrk USE input_cp2k_hfx, ONLY: create_hfx_section + USE input_cp2k_kpoints, ONLY: create_kpoint_set_section USE input_keyword_types, ONLY: keyword_create,& keyword_release,& keyword_type @@ -708,8 +709,7 @@ SUBROUTINE create_ri_g0w0(section) CALL keyword_create(keyword, __LOCATION__, name="RI_SIGMA_X", & description="If true, the exchange self-energy is calculated approximatively with RI. "// & - "This is only recommended in case exact exchange is very costly, e.g. when "// & - "using diffuse basis functions (seems not to work for periodic systems).", & + "If false, the Hartree-Fock implementation in CP2K is used.", & usage="RI_SIGMA_X", & default_l_val=.TRUE., & lone_keyword_l_val=.TRUE.) @@ -734,9 +734,11 @@ SUBROUTINE create_ri_g0w0(section) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="PERIODIC", & - description="If true, the periodic correction scheme is used employing k-points.", & - usage="PERIODIC", & + CALL keyword_create(keyword, __LOCATION__, name="PERIODIC_CORRECTION", & + description="If true, the periodic correction scheme is used employing k-points. "// & + "Method is not recommended to use, use instead PERIODIC_LOW_SCALING which much "// & + "more accurate than the periodic correction.", & + usage="PERIODIC_CORRECTION", & default_l_val=.FALSE., & lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) @@ -825,6 +827,11 @@ SUBROUTINE create_ri_g0w0(section) CALL section_add_subsection(section, subsection) CALL section_release(subsection) + ! here we generate a subsection for calculating the GW band structures + CALL create_kpoint_set_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + END SUBROUTINE create_ri_g0w0 ! ************************************************************************************************** @@ -866,8 +873,9 @@ SUBROUTINE create_periodic_gw_correction_section(section) TYPE(keyword_type), POINTER :: keyword CPASSERT(.NOT. ASSOCIATED(section)) - CALL section_create(section, __LOCATION__, name="PERIODIC", & - description="Parameters influencing correction for periodic GW.", & + CALL section_create(section, __LOCATION__, name="PERIODIC_CORRECTION", & + description="Parameters influencing correction for periodic GW. Old method, "// & + "not recommended to use", & n_keywords=12, n_subsections=1, repeats=.FALSE.) NULLIFY (keyword) @@ -1126,38 +1134,90 @@ SUBROUTINE create_low_scaling(section) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create( & + keyword, __LOCATION__, name="KPOINTS", & + description="Keyword activates periodic, low-scaling GW calculations (&LOW_SCALING section also needed). "// & + "For periodic calculations, kpoints are used for the density response, the "// & + "Coulomb interaction and the screened Coulomb interaction. For 2d periodic systems, e.g. xz "// & + "periodicity, please also specify KPOINTS, e.g. N_x 1 N_z.", & + usage="KPOINTS N_x N_y N_z", & + n_var=3, type_of_var=integer_t, default_i_vals=(/0, 0, 0/)) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + CALL keyword_create( & keyword, __LOCATION__, & - name="CUTOFF_W", & - description="Cutoff for screened Coulomb interaction for GW kpoints.", & - usage="CUTOFF_W 0.5", & - default_r_val=0.5_dp) + name="KPOINT_WEIGHTS_W", & + description="For kpoints in low-scaling GW, a Monkhorst-Pack mesh is used. The screened Coulomb "// & + "interaction W(k) needs special care near the Gamma point (e.g. in 3d, W(k) diverges at the "// & + "Gamma point with W(k) ~ k^alpha). KPOINT_WEIGHTS_W decides how the weights of the "// & + "Monkhorst-Pack mesh are chosen to compute W(R) = int_BZ W(k) exp(ikR) dk (BZ=Brllouin zone). ", & + usage="KPOINT_WEIGHTS_W AUTO", & + enum_c_vals=s2a("TAILORED", "AUTO", "UNIFORM"), & + enum_i_vals=(/kp_weights_W_tailored, kp_weights_W_auto, kp_weights_W_uniform/), & + enum_desc=s2a("Choose k-point integration weights such that the function f(k)=k^alpha is "// & + "exactly integrated. alpha is specified using EXPONENT_TAILORED_WEIGHTS.", & + "As 'TAILORED', but alpha is chosen automatically according to dimensionality "// & + "(3D: alpha = -2 for 3D, 2D: alpha = -1 for exchange self-energy, uniform "// & + "weights for correlation self-energy).", & + "Choose the same weight for every k-point (original Monkhorst-Pack method)."), & + default_i_val=kp_weights_W_auto) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) CALL keyword_create( & - keyword, __LOCATION__, name="KPOINTS", & - description="For periodic calculations, using kpoints for the density response and the "// & - "Coulomb operator are strongly recommended. For 2d periodic systems (e.g. xy "// & - "periodicity, please specify KPOINTS N_x 0 N_z.", & - usage="KPOINTS N_x N_y N_z", & - n_var=3, type_of_var=integer_t, default_i_vals=(/0, 0, 0/)) + keyword, __LOCATION__, & + name="EXPONENT_TAILORED_WEIGHTS", & + description="Gives the exponent of exactly integrated function in case 'KPOINT_WEIGHTS_W "// & + "TAILORED' is chosen.", & + usage="EXPONENT_TAILORED_WEIGHTS -2", & + default_r_val=-2.0_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create( & + keyword, __LOCATION__, name="REL_CUTOFFS_CHI_W", & + description="For periodic calculations, the irreducible density reponse chi^0(r,r') and the screened "// & + "Coulomb interaction W(r,r') is truncated "// & + "in real space. No truncation of chi^0(r,r') and W(r,r') for |r-r'| < cutoff_1*(smallest lattice vector length), "// & + "chi^0(r,r') is set to zero for |r-r'| > cutoff_2*(smallest lattice vector length). Smooth transition "// & + "in between. For cutoff_1 = 0.5, and cutoff_2 = 0.5, we have the minimum image convention (that is also default).", & + usage="REL_CUTOFFS_CHI cutoff_1 cutoff_2", & + n_var=2, type_of_var=real_t, default_r_vals=(/0.5_dp, 0.5_dp/)) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) CALL keyword_create( & keyword, __LOCATION__, & - name="EXP_KPOINTS", & - description="For kpoints in low-scaling GW, a Monkhorst-Pack mesh is used. Because the screened Coulomb "// & - "interaction W(k) diverges at the Gamma point with W(k) ~ k^alpha, we adapt the weights of the "// & - "Monkhorst-Pack mesh to compute int_BZ k^alpha dk (BZ=Brllouin zone) correctly with the Monkhorst-Pack "// & - "mesh. You can enter here the exponent alpha. For solids, the exponent is -2 (known from plane waves), "// & - "for 2d periodic systems -1 and for 1d systems W(k) ~ log(1-cos(a*k)) where a is the length of the unit "// & - "cell in periodic direction. If you enter 1.0, one of these three functions are picked according to the "// & - "periodicity. If you enter a value bigger than 2.0, the ordinary Monkhorst-Pack mesh with identical "// & - "weights is chosen.", & - usage="EXP_KPOINTS -2.0", & - default_r_val=1.0_dp) + name="REGULARIZATION_RI", & + description="Parameter to reduce the expansion coefficients in RI for periodic GW. Larger parameter "// & + "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) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create( & + keyword, __LOCATION__, & + name="MAKE_CHI_POS_DEFINITE", & + description="If true, makes eigenvalue decomposition of chi(iw,k) and removes negative "// & + "eigenvalues. May increase computational cost significantly. Only recommended to try in case "// & + "Cholesky decomposition of epsilon(iw,k) fails.", & + usage="MAKE_CHI_POS_DEFINITE", & + default_l_val=.TRUE., & + lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create( & + keyword, __LOCATION__, & + name="K_MESH_G_FACTOR", & + description="The k-mesh for the Green's function can be chosen to be larger than the k-mesh for "// & + "W (without much higher computational cost). The factor given here multiplies the mesh for W to obtain"// & + "the k-mesh for G. Example: factor 4, k-mesh for W: 4x4x1 -> k-mesh for G: 16x16x1 (z-dir. is "// & + "non-periodic).", & + default_i_val=1) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -1497,7 +1557,7 @@ SUBROUTINE create_integrals_section(section) CALL keyword_create(keyword, __LOCATION__, name="SIZE_LATTICE_SUM", & description="Size of sum range L. ", & usage="SIZE_LATTICE_SUM 10", & - default_i_val=6) + default_i_val=5) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) diff --git a/src/kpoint_coulomb_2c.F b/src/kpoint_coulomb_2c.F index 6b07f64794..90c69f6a73 100644 --- a/src/kpoint_coulomb_2c.F +++ b/src/kpoint_coulomb_2c.F @@ -61,10 +61,9 @@ MODULE kpoint_coulomb_2c !> \param atomic_kind_set ... !> \param size_lattice_sum ... !> \param operator_type ... -!> \param cutoff ... ! ************************************************************************************************** SUBROUTINE build_2c_coulomb_matrix_kp(matrix_v_kp, kpoints, basis_type, cell, particle_set, qs_kind_set, & - atomic_kind_set, size_lattice_sum, operator_type, cutoff) + atomic_kind_set, size_lattice_sum, operator_type) TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_v_kp TYPE(kpoint_type), POINTER :: kpoints @@ -74,7 +73,6 @@ SUBROUTINE build_2c_coulomb_matrix_kp(matrix_v_kp, kpoints, basis_type, cell, pa TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set INTEGER :: size_lattice_sum, operator_type - REAL(dp) :: cutoff CHARACTER(LEN=*), PARAMETER :: routineN = 'build_2c_coulomb_matrix_kp' @@ -89,7 +87,7 @@ SUBROUTINE build_2c_coulomb_matrix_kp(matrix_v_kp, kpoints, basis_type, cell, pa CALL lattice_sum(matrix_v_kp, kpoints, basis_type, cell, particle_set, & qs_kind_set, atomic_kind_set, size_lattice_sum, matrix_v_L_tmp, & - operator_type, cutoff) + operator_type) CALL deallocate_tmp(matrix_v_L_tmp) @@ -109,11 +107,10 @@ END SUBROUTINE build_2c_coulomb_matrix_kp !> \param size_lattice_sum ... !> \param matrix_v_L_tmp ... !> \param operator_type ... -!> \param cutoff ... ! ************************************************************************************************** SUBROUTINE lattice_sum(matrix_v_kp, kpoints, basis_type, cell, particle_set, & qs_kind_set, atomic_kind_set, size_lattice_sum, matrix_v_L_tmp, & - operator_type, cutoff) + operator_type) TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_v_kp TYPE(kpoint_type), POINTER :: kpoints @@ -125,7 +122,6 @@ SUBROUTINE lattice_sum(matrix_v_kp, kpoints, basis_type, cell, particle_set, & INTEGER :: size_lattice_sum TYPE(dbcsr_type), POINTER :: matrix_v_L_tmp INTEGER :: operator_type - REAL(dp) :: cutoff CHARACTER(LEN=*), PARAMETER :: routineN = 'lattice_sum' @@ -213,7 +209,7 @@ SUBROUTINE lattice_sum(matrix_v_kp, kpoints, basis_type, cell, particle_set, & ! Compute (P 0 | Q vec_L) and store it in matrix_v_L_tmp CALL compute_v_transl(matrix_v_L_tmp, blocks_of_v_L, vec_L, particle_set, & qs_kind_set, atomic_kind_set, basis_type, cell, & - operator_type, cutoff) + operator_type) DO i_block = 1, SIZE(blocks_of_v_L) blocks_of_v_L_store(i_block)%block(:, :) = blocks_of_v_L_store(i_block)%block(:, :) & @@ -320,10 +316,9 @@ SUBROUTINE set_blocks_to_matrix_v_kp(matrix_v_kp, blocks_of_v_kp) !> \param basis_type ... !> \param cell ... !> \param operator_type ... -!> \param cutoff ... ! ************************************************************************************************** SUBROUTINE compute_v_transl(matrix_v_L_tmp, blocks_of_v_L, vec_L, particle_set, & - qs_kind_set, atomic_kind_set, basis_type, cell, operator_type, cutoff) + qs_kind_set, atomic_kind_set, basis_type, cell, operator_type) TYPE(dbcsr_type), POINTER :: matrix_v_L_tmp TYPE(two_d_util_type), ALLOCATABLE, DIMENSION(:) :: blocks_of_v_L @@ -334,7 +329,6 @@ SUBROUTINE compute_v_transl(matrix_v_L_tmp, blocks_of_v_L, vec_L, particle_set, CHARACTER(LEN=*), INTENT(IN) :: basis_type TYPE(cell_type), POINTER :: cell INTEGER :: operator_type - REAL(dp) :: cutoff CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_v_transl' @@ -384,7 +378,7 @@ SUBROUTINE compute_v_transl(matrix_v_L_tmp, blocks_of_v_L, vec_L, particle_set, blocks_of_v_L(i_block)%block = 0.0_dp CALL int_operators_r12_ab_shg(operator_type, blocks_of_v_L(i_block)%block, rab=rab_L, & - fba=basis_set_a, fbb=basis_set_b, scona_shg=contr_a, sconb_shg=contr_b, omega=cutoff, & + fba=basis_set_a, fbb=basis_set_b, scona_shg=contr_a, sconb_shg=contr_b, & calculate_forces=.FALSE.) i_block = i_block + 1 diff --git a/src/kpoint_methods.F b/src/kpoint_methods.F index bede6f58d6..ec4c10e06b 100644 --- a/src/kpoint_methods.F +++ b/src/kpoint_methods.F @@ -73,6 +73,7 @@ MODULE kpoint_methods qs_matrix_pools_type USE qs_mo_types, ONLY: allocate_mo_set,& get_mo_set,& + init_mo_set,& mo_set_p_type,& mo_set_type,& set_mo_set @@ -93,7 +94,7 @@ MODULE kpoint_methods CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpoint_methods' - PUBLIC :: kpoint_initialize, kpoint_env_initialize, kpoint_initialize_mos + PUBLIC :: kpoint_initialize, kpoint_env_initialize, kpoint_initialize_mos, kpoint_initialize_mo_set PUBLIC :: kpoint_init_cell_index, kpoint_set_mo_occupation PUBLIC :: kpoint_density_matrices, kpoint_density_transform PUBLIC :: rskp_transform @@ -548,6 +549,42 @@ SUBROUTINE kpoint_initialize_mos(kpoint, mos, added_mos) END SUBROUTINE kpoint_initialize_mos +! ************************************************************************************************** +!> \brief ... +!> \param kpoint ... +! ************************************************************************************************** + SUBROUTINE kpoint_initialize_mo_set(kpoint) + TYPE(kpoint_type), POINTER :: kpoint + + CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_initialize_mo_set' + + INTEGER :: handle, ic, ik, ikk, ispin + TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools + TYPE(cp_fm_type), POINTER :: mo_coeff + TYPE(mo_set_p_type), DIMENSION(:, :), POINTER :: moskp + + CALL timeset(routineN, handle) + + DO ik = 1, SIZE(kpoint%kp_env) + CALL mpools_get(kpoint%mpools, ao_mo_fm_pools=ao_mo_fm_pools) + moskp => kpoint%kp_env(ik)%kpoint_env%mos + ikk = kpoint%kp_range(1) + ik - 1 + CPASSERT(ASSOCIATED(moskp)) + DO ispin = 1, SIZE(moskp, 2) + DO ic = 1, SIZE(moskp, 1) + CALL get_mo_set(moskp(ic, ispin)%mo_set, mo_coeff=mo_coeff) + IF (.NOT. ASSOCIATED(mo_coeff)) THEN + CALL init_mo_set(moskp(ic, ispin)%mo_set, & + fm_pool=ao_mo_fm_pools(ispin)%pool, name="kpoints") + END IF + END DO + END DO + END DO + + CALL timestop(handle) + + END SUBROUTINE kpoint_initialize_mo_set + ! ************************************************************************************************** !> \brief Generates the mapping of cell indices and linear RS index !> CELL (0,0,0) is always mapped to index 1 diff --git a/src/mp2_integrals.F b/src/mp2_integrals.F index 733136032d..2bdd354ffb 100644 --- a/src/mp2_integrals.F +++ b/src/mp2_integrals.F @@ -53,14 +53,10 @@ MODULE mp2_integrals USE hfx_types, ONLY: alloc_containers,& block_ind_type,& hfx_compression_type - USE input_constants, ONLY: do_eri_gpw,& - do_eri_mme,& - do_eri_os,& - do_potential_coulomb,& - do_potential_id,& - do_potential_long,& - do_potential_short,& - do_potential_truncated + USE input_constants, ONLY: & + do_eri_gpw, do_eri_mme, do_eri_os, do_potential_coulomb, do_potential_id, & + do_potential_long, do_potential_short, do_potential_truncated, kp_weights_W_auto, & + kp_weights_W_tailored, kp_weights_W_uniform USE kinds, ONLY: default_string_length,& dp,& int_8 @@ -1252,13 +1248,11 @@ SUBROUTINE compute_kpoints(qs_env, kpoints, unit_nr) kpoints%eps_geo = 1.e-6_dp kpoints%full_grid = .TRUE. nkp_grid(1:3) = qs_env%mp2_env%ri_rpa_im_time%kp_grid(1:3) - kpoints%nkp_grid(1:3) = nkp_grid(1:3) - num_dim = periodic(1) + periodic(2) + periodic(3) DO i_dim = 1, 3 IF (periodic(i_dim) == 1) THEN - CPASSERT(MODULO(kpoints%nkp_grid(i_dim), 2) == 0) + CPASSERT(MODULO(nkp_grid(i_dim), 2) == 0) END IF END DO @@ -1270,12 +1264,10 @@ SUBROUTINE compute_kpoints(qs_env, kpoints, unit_nr) nkp = MAX(nkp_grid(1)/2, nkp_grid(2)/2, nkp_grid(3)/2) END IF - kpoints%nkp = nkp - - IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & - "KPOINT_INFO| Number of kpoints for V and W:", nkp - ALLOCATE (kpoints%xkp(3, nkp), kpoints%wkp(nkp)) + + kpoints%nkp_grid(1:3) = nkp_grid(1:3) + kpoints%nkp = nkp kpoints%wkp(:) = 1.0_dp/REAL(nkp, KIND=dp) i = 0 @@ -1298,6 +1290,25 @@ SUBROUTINE compute_kpoints(qs_env, kpoints, unit_nr) CALL set_qs_env(qs_env, kpoints=kpoints) + IF (unit_nr > 0) THEN + + WRITE (UNIT=unit_nr, FMT="(T3,A,T69,3I4)") "KPOINT_INFO| K-point mesh for V and W:", nkp_grid(1:3) + WRITE (UNIT=unit_nr, FMT="(T3,A,T75,I6)") "KPOINT_INFO| Number of kpoints for V and W:", nkp + + SELECT CASE (qs_env%mp2_env%ri_rpa_im_time%kpoint_weights_W_method) + CASE (kp_weights_W_tailored) + WRITE (UNIT=unit_nr, FMT="(T3,A,T81)") & + "KPOINT_INFO| K-point weights for W: TAILORED" + CASE (kp_weights_W_auto) + WRITE (UNIT=unit_nr, FMT="(T3,A,T81)") & + "KPOINT_INFO| K-point weights for W: AUTO" + CASE (kp_weights_W_uniform) + WRITE (UNIT=unit_nr, FMT="(T3,A,T81)") & + "KPOINT_INFO| K-point weights for W: UNIFORM" + END SELECT + + END IF + CALL timestop(handle) END SUBROUTINE compute_kpoints diff --git a/src/mp2_ri_2c.F b/src/mp2_ri_2c.F index b05e615472..c9203bd7df 100644 --- a/src/mp2_ri_2c.F +++ b/src/mp2_ri_2c.F @@ -22,12 +22,9 @@ MODULE mp2_ri_2c USE cp_blacs_env, ONLY: cp_blacs_env_create,& cp_blacs_env_release,& cp_blacs_env_type - USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_decompose,& - cp_cfm_gemm,& - cp_cfm_scale_and_add_fm,& - cp_cfm_triangular_invert + USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm,& + cp_cfm_scale_and_add_fm USE cp_cfm_types, ONLY: cp_cfm_create,& - cp_cfm_get_info,& cp_cfm_release,& cp_cfm_to_fm,& cp_cfm_type @@ -97,6 +94,7 @@ MODULE mp2_ri_2c USE qs_tensors, ONLY: build_2c_integrals,& build_2c_neighbor_lists USE rpa_communication, ONLY: communicate_buffer + USE rpa_kpoints, ONLY: cp_cfm_robust_cholesky !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num #include "./base/base_uses.f90" @@ -242,6 +240,8 @@ SUBROUTINE get_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_su CALL compute_V_by_lattice_sum(qs_env, fm_matrix_L_kpoints, fm_matrix_L_RI_metric, kpoints) + CALL regularize_RI_expansion(fm_matrix_L_RI_metric, qs_env%mp2_env%ri_rpa_im_time%regularization_RI) + CALL inversion_of_S_and_mult_with_chol_dec_of_V(fm_matrix_L_RI_metric, fm_matrix_L_kpoints, & dimen_RI, kpoints) @@ -382,8 +382,6 @@ SUBROUTINE compute_V_by_lattice_sum(qs_env, fm_matrix_L_kpoints, fm_matrix_L_RI_ INTEGER :: handle, i_dim, i_real_imag, ikp, nkp INTEGER, DIMENSION(3) :: periodic - REAL(dp) :: cutoff - REAL(KIND=dp), DIMENSION(3) :: abc TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cell_type), POINTER :: cell TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_RI_aux_transl, matrix_v_RI_kp @@ -440,23 +438,11 @@ SUBROUTINE compute_V_by_lattice_sum(qs_env, fm_matrix_L_kpoints, fm_matrix_L_RI_ END DO - CALL get_cell(cell=cell, abc=abc, periodic=periodic) - IF (periodic(1) == 0) THEN - abc(1) = 1.0E10_dp - END IF - IF (periodic(2) == 0) THEN - abc(2) = 1.0E10_dp - END IF - IF (periodic(3) == 0) THEN - abc(3) = 1.0E10_dp - END IF - cutoff = 1.0_dp/(qs_env%mp2_env%ri_rpa_im_time%cutoff*MIN(abc(1), abc(2), abc(3))) - CALL build_2c_coulomb_matrix_kp(matrix_v_RI_kp, kpoints, basis_type="RI_AUX", & cell=cell, particle_set=particle_set, qs_kind_set=qs_kind_set, & atomic_kind_set=atomic_kind_set, & size_lattice_sum=qs_env%mp2_env%mp2_gpw%size_lattice_sum, & - operator_type=operator_coulomb, cutoff=cutoff) + operator_type=operator_coulomb) DO ikp = 1, nkp @@ -1379,6 +1365,50 @@ SUBROUTINE cholesky_decomp(fm_matrix_L, dimen_RI, do_inversion) END SUBROUTINE cholesky_decomp +! ************************************************************************************************** +!> \brief ... +!> \param fm_matrix_L_RI_metric ... +!> \param regularization_RI ... +! ************************************************************************************************** + SUBROUTINE regularize_RI_expansion(fm_matrix_L_RI_metric, regularization_RI) + TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_matrix_L_RI_metric + REAL(KIND=dp) :: regularization_RI + + CHARACTER(LEN=*), PARAMETER :: routineN = 'regularize_RI_expansion' + + INTEGER :: handle, i_global, iiB, ikp, j_global, & + jjB, ncol_local, nkp, nrow_local + INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices + + CALL timeset(routineN, handle) + + nkp = SIZE(fm_matrix_L_RI_metric, 1) + + CALL cp_fm_get_info(matrix=fm_matrix_L_RI_metric(1, 1)%matrix, & + nrow_local=nrow_local, & + ncol_local=ncol_local, & + row_indices=row_indices, & + col_indices=col_indices) + + DO ikp = 1, nkp + + DO iiB = 1, nrow_local + i_global = row_indices(iiB) + DO jjB = 1, ncol_local + j_global = col_indices(jjB) + IF (j_global == i_global) THEN + fm_matrix_L_RI_metric(ikp, 1)%matrix%local_data(iiB, jjB) = & + fm_matrix_L_RI_metric(ikp, 1)%matrix%local_data(iiB, jjB) + regularization_RI + END IF + END DO + END DO + + END DO + + CALL timestop(handle) + + END SUBROUTINE regularize_RI_expansion + ! ************************************************************************************************** !> \brief ... !> \param fm_matrix_L_RI_metric ... @@ -1398,10 +1428,7 @@ SUBROUTINE inversion_of_S_and_mult_with_chol_dec_of_V(fm_matrix_L_RI_metric, fm_ 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, i_global, iiB, ikp, info_chol, & - j_global, jjB, ncol_local, nkp, & - nrow_local - INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices + INTEGER :: handle, ikp, nkp TYPE(cp_cfm_type), POINTER :: cfm_matrix_K_tmp, cfm_matrix_L_tmp, & cfm_matrix_S_inv_tmp, cfm_matrix_S_tmp TYPE(cp_fm_struct_type), POINTER :: matrix_struct @@ -1427,49 +1454,13 @@ SUBROUTINE inversion_of_S_and_mult_with_chol_dec_of_V(fm_matrix_L_RI_metric, fm_ CALL cp_cfm_scale_and_add_fm(czero, cfm_matrix_L_tmp, cone, fm_matrix_L_kpoints(ikp, 1)%matrix) CALL cp_cfm_scale_and_add_fm(cone, cfm_matrix_L_tmp, ione, fm_matrix_L_kpoints(ikp, 2)%matrix) - ! Cholesky decomposition S=U(T)U - CALL cp_cfm_cholesky_decompose(cfm_matrix_S_tmp, info_out=info_chol) - CPASSERT(info_chol == 0) - - ! Invert to get U^(-1) - CALL cp_cfm_triangular_invert(cfm_matrix_S_tmp, uplo='U') + CALL cp_cfm_robust_cholesky(cfm_matrix_S_tmp, cfm_matrix_S_inv_tmp, threshold=1.0E-12_dp, exponent=-0.5_dp) - ! clean the lower part of the U^{-1} matrix (just to not have surprises afterwards) - CALL cp_cfm_get_info(matrix=cfm_matrix_S_tmp, & - nrow_local=nrow_local, & - ncol_local=ncol_local, & - row_indices=row_indices, & - col_indices=col_indices) - - DO iiB = 1, nrow_local - i_global = row_indices(iiB) - DO jjB = 1, ncol_local - j_global = col_indices(jjB) - IF (j_global < i_global) cfm_matrix_S_tmp%local_data(iiB, jjB) = czero - END DO - END DO - - ! get S^(-1) = U^(-1)U^(-H) - CALL cp_cfm_gemm("N", "C", dimen_RI, dimen_RI, dimen_RI, cone, cfm_matrix_S_tmp, cfm_matrix_S_tmp, & + ! 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, & czero, cfm_matrix_S_inv_tmp) - ! Cholesky decomposition V=L(T)L - CALL cp_cfm_cholesky_decompose(cfm_matrix_L_tmp, info_out=info_chol) - CPASSERT(info_chol == 0) - - ! clean the lower part of the L matrix (just to not have surprises afterwards) - CALL cp_cfm_get_info(matrix=cfm_matrix_L_tmp, & - nrow_local=nrow_local, & - ncol_local=ncol_local, & - row_indices=row_indices, & - col_indices=col_indices) - DO iiB = 1, nrow_local - i_global = row_indices(iiB) - DO jjB = 1, ncol_local - j_global = col_indices(jjB) - IF (j_global < i_global) cfm_matrix_L_tmp%local_data(iiB, jjB) = czero - END DO - END DO + CALL cp_cfm_robust_cholesky(cfm_matrix_L_tmp, cfm_matrix_S_tmp, threshold=0.0_dp, exponent=0.5_dp) ! get K = S^(-1)*L, where L is the Cholesky decomposition of V = L^T*L CALL cp_cfm_gemm("N", "C", dimen_RI, dimen_RI, dimen_RI, cone, cfm_matrix_S_inv_tmp, cfm_matrix_L_tmp, & @@ -1568,4 +1559,5 @@ SUBROUTINE create_matrix_L(fm_matrix_L, blacs_env_L, dimen_RI, para_env_L, name, CALL timestop(handle) END SUBROUTINE create_matrix_L + END MODULE mp2_ri_2c diff --git a/src/mp2_setup.F b/src/mp2_setup.F index cbd7f660e9..5b76d3e4f0 100644 --- a/src/mp2_setup.F +++ b/src/mp2_setup.F @@ -17,13 +17,16 @@ MODULE mp2_setup cp_logger_type USE cp_output_handling, ONLY: cp_print_key_finished_output,& cp_print_key_unit_nr + USE cp_parser_methods, ONLY: read_float_object USE input_constants, ONLY: & do_eri_mme, do_potential_short, mp2_method_direct, mp2_method_gpw, mp2_method_none, & mp2_ri_optimize_basis, ri_mp2_laplace, ri_mp2_method_gpw, ri_rpa_method_gpw USE input_section_types, ONLY: section_vals_get_subs_vals,& section_vals_type,& section_vals_val_get - USE kinds, ONLY: dp + USE kinds, ONLY: default_string_length,& + dp,& + max_line_length USE machine, ONLY: m_flush USE mathlib, ONLY: erfc_cutoff USE mp2_types, ONLY: mp2_method_direct,& @@ -57,7 +60,10 @@ SUBROUTINE read_mp2_section(input, mp2_env) CHARACTER(len=*), PARAMETER :: routineN = 'read_mp2_section' - INTEGER :: handle, ival, unit_nr + CHARACTER(LEN=default_string_length), & + DIMENSION(:), POINTER :: string_pointer + CHARACTER(LEN=max_line_length) :: error_message + INTEGER :: handle, i, i_special_kp, ival, unit_nr INTEGER, DIMENSION(:), POINTER :: tmplist LOGICAL :: do_mp2, do_opt_ri_basis, do_ri_mp2, & do_ri_sos_mp2, do_rpa @@ -163,29 +169,46 @@ SUBROUTINE read_mp2_section(input, mp2_env) CALL section_vals_val_get(mp2_section, "RI_RPA%GW%IC%EPS_DIST", & r_val=mp2_env%ri_g0w0%eps_dist) - CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC", & + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC_CORRECTION", & l_val=mp2_env%ri_g0w0%do_periodic) - CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC%KPOINTS", & + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC_CORRECTION%KPOINTS", & i_vals=mp2_env%ri_g0w0%kp_grid) - CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC%NUM_KP_GRIDS", & + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC_CORRECTION%NUM_KP_GRIDS", & i_val=mp2_env%ri_g0w0%num_kp_grids) - CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC%EPS_KPOINT", & + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC_CORRECTION%EPS_KPOINT", & r_val=mp2_env%ri_g0w0%eps_kpoint) - CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC%MO_COEFF_GAMMA", & + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC_CORRECTION%MO_COEFF_GAMMA", & l_val=mp2_env%ri_g0w0%do_mo_coeff_gamma) - CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC%AVERAGE_DEGENERATE_LEVELS", & + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC_CORRECTION%AVERAGE_DEGENERATE_LEVELS", & l_val=mp2_env%ri_g0w0%do_average_deg_levels) - CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC%EPS_EIGENVAL", & + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC_CORRECTION%EPS_EIGENVAL", & r_val=mp2_env%ri_g0w0%eps_eigenval) - CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC%EXTRAPOLATE_KPOINTS", & + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC_CORRECTION%EXTRAPOLATE_KPOINTS", & l_val=mp2_env%ri_g0w0%do_extra_kpoints) - CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC%DO_AUX_BAS_GW", & + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC_CORRECTION%DO_AUX_BAS_GW", & l_val=mp2_env%ri_g0w0%do_aux_bas_gw) - CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC%FRACTION_AUX_MOS", & + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC_CORRECTION%FRACTION_AUX_MOS", & r_val=mp2_env%ri_g0w0%frac_aux_mos) - CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC%NUM_OMEGA_POINTS", & + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PERIODIC_CORRECTION%NUM_OMEGA_POINTS", & i_val=mp2_env%ri_g0w0%num_omega_points) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%KPOINT_SET%NPOINTS", & + i_val=mp2_env%ri_g0w0%n_kp_in_kp_line) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%KPOINT_SET%SPECIAL_POINT", & + n_rep_val=mp2_env%ri_g0w0%n_special_kp) + ALLOCATE (mp2_env%ri_g0w0%xkp_special_kp(3, mp2_env%ri_g0w0%n_special_kp)) + DO i_special_kp = 1, mp2_env%ri_g0w0%n_special_kp + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%KPOINT_SET%SPECIAL_POINT", & + i_rep_val=i_special_kp, c_vals=string_pointer) + CPASSERT(SIZE(string_pointer(:), 1) == 4) + DO i = 1, 3 + CALL read_float_object(string_pointer(i + 1), & + mp2_env%ri_g0w0%xkp_special_kp(i, i_special_kp), & + error_message) + IF (LEN_TRIM(error_message) > 0) CPABORT(TRIM(error_message)) + END DO + END DO + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%PRINT%SELF_ENERGY", & l_val=mp2_env%ri_g0w0%print_self_energy) @@ -204,12 +227,20 @@ SUBROUTINE read_mp2_section(input, mp2_env) CALL section_vals_val_get(low_scaling_section, "DO_KPOINTS", & l_val=mp2_env%ri_rpa_im_time%do_im_time_kpoints) - CALL section_vals_val_get(low_scaling_section, "CUTOFF_W", & - r_val=mp2_env%ri_rpa_im_time%cutoff) CALL section_vals_val_get(low_scaling_section, "KPOINTS", & i_vals=mp2_env%ri_rpa_im_time%kp_grid) - CALL section_vals_val_get(low_scaling_section, "EXP_KPOINTS", & - r_val=mp2_env%ri_rpa_im_time%exp_kpoints) + 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", & + r_val=mp2_env%ri_rpa_im_time%exp_tailored_weights) + CALL section_vals_val_get(low_scaling_section, "REL_CUTOFFS_CHI_W", & + 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, "MAKE_CHI_POS_DEFINITE", & + l_val=mp2_env%ri_rpa_im_time%make_chi_pos_definite) + 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, "MIN_BLOCK_SIZE", & i_val=mp2_env%ri_rpa_im_time%min_bsize) diff --git a/src/mp2_types.F b/src/mp2_types.F index 3e483b9edd..f9e686a7c7 100644 --- a/src/mp2_types.F +++ b/src/mp2_types.F @@ -30,6 +30,7 @@ MODULE mp2_types ri_mp2_method_gpw,& ri_rpa_method_gpw USE kinds, ONLY: dp + USE kpoint_types, ONLY: kpoint_type USE libint_2c_3c, ONLY: libint_potential_type USE machine, ONLY: m_walltime USE message_passing, ONLY: mp_sum @@ -143,13 +144,16 @@ MODULE mp2_types TYPE ri_rpa_im_time_type INTEGER :: cut_memory - LOGICAL :: memory_info - REAL(KIND=dp) :: eps_filter, cutoff, & - exp_kpoints, eps_filter_factor, eps_compress - INTEGER :: group_size_P, group_size_3c + LOGICAL :: memory_info, make_chi_pos_definite + REAL(KIND=dp) :: eps_filter, & + eps_filter_factor, eps_compress, exp_tailored_weights, regularization_RI + REAL(KIND=dp), DIMENSION(:), POINTER :: rel_cutoffs_chi_W + REAL(KIND=dp), DIMENSION(2) :: abs_cutoffs_chi_W + INTEGER :: group_size_P, group_size_3c, kpoint_weights_W_method, k_mesh_g_factor INTEGER, DIMENSION(:), POINTER :: kp_grid LOGICAL :: do_im_time_kpoints INTEGER :: min_bsize, min_bsize_mo + TYPE(kpoint_type), POINTER :: kpoints_G, kpoints_Sigma, kpoints_Sigma_no_xc END TYPE TYPE ri_g0w0_type @@ -186,6 +190,8 @@ MODULE mp2_types INTEGER :: print_exx LOGICAL :: do_gamma_only_sigma LOGICAL :: update_xc_energy + INTEGER :: n_kp_in_kp_line, n_special_kp + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: xkp_special_kp END TYPE TYPE ri_basis_opt diff --git a/src/qs_band_structure.F b/src/qs_band_structure.F index dcd7dad339..5ce07120b0 100644 --- a/src/qs_band_structure.F +++ b/src/qs_band_structure.F @@ -12,15 +12,12 @@ !> \author JGH ! ************************************************************************************************** MODULE qs_band_structure - USE cell_types, ONLY: cell_type,& - get_cell + USE cell_types, ONLY: cell_type USE cp_blacs_env, ONLY: cp_blacs_env_retain,& cp_blacs_env_type USE cp_control_types, ONLY: dft_control_type USE cp_files, ONLY: close_file,& open_file - USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type - USE cp_fm_types, ONLY: cp_fm_type USE cp_log_handling, ONLY: cp_logger_get_default_io_unit USE cp_para_env, ONLY: cp_para_env_retain USE cp_para_types, ONLY: cp_para_env_type @@ -35,6 +32,7 @@ MODULE qs_band_structure max_line_length USE kpoint_methods, ONLY: kpoint_env_initialize,& kpoint_init_cell_index,& + kpoint_initialize_mo_set,& kpoint_initialize_mos USE kpoint_types, ONLY: get_kpoint_info,& kpoint_create,& @@ -51,10 +49,7 @@ MODULE qs_band_structure qs_env_release,& qs_environment_type USE qs_gamma2kp, ONLY: create_kp_from_gamma - USE qs_matrix_pools, ONLY: mpools_get - USE qs_mo_types, ONLY: get_mo_set,& - init_mo_set,& - mo_set_p_type + USE qs_mo_types, ONLY: get_mo_set USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type USE qs_scf_diagonalization, ONLY: do_general_diag_kp USE qs_scf_types, ONLY: qs_scf_env_type @@ -315,18 +310,11 @@ SUBROUTINE calculate_kp_orbitals(qs_env, kpoint, scheme, nadd, mp_grid, kpgenera OPTIONAL :: kpgeneral INTEGER, INTENT(IN), OPTIONAL :: group_size_ext - INTEGER :: i, ic, ik, ikk, ispin, ix, iy, iz, & - npoints - REAL(KIND=dp), DIMENSION(3) :: kpt_latt - REAL(KIND=dp), DIMENSION(3, 3) :: h_inv - TYPE(cell_type), POINTER :: cell + INTEGER :: i, ix, iy, iz, npoints TYPE(cp_blacs_env_type), POINTER :: blacs_env - TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools - TYPE(cp_fm_type), POINTER :: mo_coeff TYPE(cp_para_env_type), POINTER :: para_env TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s TYPE(dft_control_type), POINTER :: dft_control - TYPE(mo_set_p_type), DIMENSION(:, :), POINTER :: mos TYPE(neighbor_list_set_p_type), DIMENSION(:), & POINTER :: sab_nl TYPE(qs_scf_env_type), POINTER :: scf_env @@ -365,17 +353,14 @@ SUBROUTINE calculate_kp_orbitals(qs_env, kpoint, scheme, nadd, mp_grid, kpgenera kpoint%nkp = npoints ALLOCATE (kpoint%xkp(3, npoints), kpoint%wkp(npoints)) kpoint%wkp(:) = 1._dp/REAL(npoints, KIND=dp) - CALL get_qs_env(qs_env, cell=cell) - CALL get_cell(cell, h_inv=h_inv) i = 0 DO ix = 1, mp_grid(1) DO iy = 1, mp_grid(2) DO iz = 1, mp_grid(3) i = i + 1 - kpt_latt(1) = REAL(2*ix - mp_grid(1) - 1, KIND=dp)/(2._dp*REAL(mp_grid(1), KIND=dp)) - kpt_latt(2) = REAL(2*iy - mp_grid(2) - 1, KIND=dp)/(2._dp*REAL(mp_grid(2), KIND=dp)) - kpt_latt(3) = REAL(2*iz - mp_grid(3) - 1, KIND=dp)/(2._dp*REAL(mp_grid(3), KIND=dp)) - kpoint%xkp(1:3, i) = MATMUL(TRANSPOSE(h_inv), kpt_latt(:)) + kpoint%xkp(1, i) = REAL(2*ix - mp_grid(1) - 1, KIND=dp)/(2._dp*REAL(mp_grid(1), KIND=dp)) + kpoint%xkp(2, i) = REAL(2*iy - mp_grid(2) - 1, KIND=dp)/(2._dp*REAL(mp_grid(2), KIND=dp)) + kpoint%xkp(3, i) = REAL(2*iz - mp_grid(3) - 1, KIND=dp)/(2._dp*REAL(mp_grid(3), KIND=dp)) END DO END DO END DO @@ -412,24 +397,11 @@ SUBROUTINE calculate_kp_orbitals(qs_env, kpoint, scheme, nadd, mp_grid, kpgenera CALL kpoint_env_initialize(kpoint) CALL kpoint_initialize_mos(kpoint, qs_env%mos, nadd) - DO ik = 1, SIZE(kpoint%kp_env) - CALL mpools_get(kpoint%mpools, ao_mo_fm_pools=ao_mo_fm_pools) - mos => kpoint%kp_env(ik)%kpoint_env%mos - ikk = kpoint%kp_range(1) + ik - 1 - CPASSERT(ASSOCIATED(mos)) - DO ispin = 1, SIZE(mos, 2) - DO ic = 1, SIZE(mos, 1) - CALL get_mo_set(mos(ic, ispin)%mo_set, mo_coeff=mo_coeff) - IF (.NOT. ASSOCIATED(mo_coeff)) THEN - CALL init_mo_set(mos(ic, ispin)%mo_set, & - fm_pool=ao_mo_fm_pools(ispin)%pool, name="kpoints") - END IF - END DO - END DO - END DO + 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.) diff --git a/src/qs_gamma2kp.F b/src/qs_gamma2kp.F index 1507c12747..80b85f7ff6 100644 --- a/src/qs_gamma2kp.F +++ b/src/qs_gamma2kp.F @@ -14,10 +14,18 @@ MODULE qs_gamma2kp USE cp_para_types, ONLY: cp_para_env_type USE cp_subsys_types, ONLY: cp_subsys_type - USE input_section_types, ONLY: section_vals_type + USE dbcsr_api, ONLY: dbcsr_add,& + dbcsr_p_type + USE input_constants, ONLY: atomic_guess,& + xc_none + USE input_section_types, ONLY: section_vals_get_subs_vals,& + section_vals_type,& + section_vals_val_get,& + section_vals_val_set USE kinds, ONLY: dp USE kpoint_types, ONLY: kpoint_create,& kpoint_type + USE pw_types, ONLY: pw_p_type USE qs_energy_init, ONLY: qs_energies_init USE qs_environment, ONLY: qs_init USE qs_environment_types, ONLY: get_qs_env,& @@ -25,9 +33,12 @@ MODULE qs_gamma2kp qs_environment_type,& set_qs_env USE qs_ks_methods, ONLY: qs_ks_update_qs_env + USE qs_rho_types, ONLY: qs_rho_get,& + qs_rho_type USE qs_scf_initialization, ONLY: qs_scf_env_init_basic USE qs_scf_types, ONLY: qs_scf_env_type,& scf_env_release + USE scf_control_types, ONLY: scf_control_type #include "./base/base_uses.f90" IMPLICIT NONE @@ -45,15 +56,25 @@ MODULE qs_gamma2kp !> \brief ... !> \param qs_env ... !> \param qs_env_kp ... +!> \param with_xc_terms ... ! ************************************************************************************************** - SUBROUTINE create_kp_from_gamma(qs_env, qs_env_kp) + SUBROUTINE create_kp_from_gamma(qs_env, qs_env_kp, with_xc_terms) TYPE(qs_environment_type), POINTER :: qs_env, qs_env_kp + LOGICAL, OPTIONAL :: with_xc_terms + INTEGER :: ispin, xc_func + LOGICAL :: without_xc_terms TYPE(cp_para_env_type), POINTER :: para_env TYPE(cp_subsys_type), POINTER :: cp_subsys + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp, rho_ao_kp_gamma TYPE(kpoint_type), POINTER :: kpoint + TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_g_gamma, rho_g_kp, rho_r_gamma, & + rho_r_kp + TYPE(qs_rho_type), POINTER :: rho, rho_gamma TYPE(qs_scf_env_type), POINTER :: scf_env - TYPE(section_vals_type), POINTER :: force_env_section, subsys_section + TYPE(scf_control_type), POINTER :: scf_control + TYPE(section_vals_type), POINTER :: force_env_section, subsys_section, & + xc_fun_section !--------------------------------------------------------------------------------------------! @@ -76,20 +97,50 @@ SUBROUTINE create_kp_from_gamma(qs_env, qs_env_kp) kpoint%use_real_wfn = .TRUE. kpoint%parallel_group_size = 0 + without_xc_terms = .FALSE. + IF (PRESENT(with_xc_terms)) without_xc_terms = .NOT. with_xc_terms + + IF (without_xc_terms) THEN + xc_fun_section => section_vals_get_subs_vals(force_env_section, "DFT%XC%XC_FUNCTIONAL") + CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=xc_func) + CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", i_val=xc_none) + END IF + CALL qs_env_create(qs_env_kp) CALL qs_init(qs_env_kp, para_env, cp_subsys=cp_subsys, kpoint_env=kpoint, & force_env_section=force_env_section, subsys_section=subsys_section, & use_motion_section=.FALSE.) + CALL get_qs_env(qs_env_kp, scf_control=scf_control) + scf_control%density_guess = atomic_guess CALL qs_energies_init(qs_env_kp, calc_forces=.FALSE.) NULLIFY (scf_env) CALL qs_scf_env_init_basic(qs_env_kp, scf_env) + CALL set_qs_env(qs_env_kp, scf_env=scf_env) CALL scf_env_release(scf_env) + ! copy density matrix, n(r) and n(G) from Gamma-only to kpoint calculation + CALL get_qs_env(qs_env, rho=rho_gamma) + CALL qs_rho_get(rho_gamma, rho_ao_kp=rho_ao_kp_gamma, rho_r=rho_r_gamma, rho_g=rho_g_gamma) + CALL get_qs_env(qs_env_kp, rho=rho) + CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp, rho_r=rho_r_kp, rho_g=rho_g_kp) + + DO ispin = 1, SIZE(rho_r_gamma) + rho_r_kp(ispin)%pw%cr3d(:, :, :) = rho_r_gamma(ispin)%pw%cr3d(:, :, :) + rho_g_kp(ispin)%pw%cc(:) = rho_g_gamma(ispin)%pw%cc(:) + CALL dbcsr_add(matrix_a=rho_ao_kp(ispin, 1)%matrix, alpha_scalar=0.0_dp, & + matrix_b=rho_ao_kp_gamma(ispin, 1)%matrix, beta_scalar=1.0_dp) + END DO + CALL qs_ks_update_qs_env(qs_env_kp, print_active=.FALSE.) + IF (without_xc_terms) THEN + ! set back the functional + CALL section_vals_val_set(xc_fun_section, "_SECTION_PARAMETERS_", i_val=xc_func) + END IF + END SUBROUTINE create_kp_from_gamma END MODULE qs_gamma2kp diff --git a/src/qs_ks_methods.F b/src/qs_ks_methods.F index 6fe464b889..890a7becce 100644 --- a/src/qs_ks_methods.F +++ b/src/qs_ks_methods.F @@ -678,6 +678,7 @@ SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, & ! sum up potentials and integrate ! Pointing my_rho to the density matrix rho_ao my_rho => rho_ao + CALL sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, vppl_rspace, & v_rspace_new, v_rspace_new_aux_fit, v_tau_rspace, v_tau_rspace_aux_fit, & v_efield_rspace, v_sic_rspace, v_spin_ddapc_rest_r, v_sccs_rspace, v_rspace_embed, & diff --git a/src/qs_wannier90.F b/src/qs_wannier90.F index 4204c3a6b7..ffad942192 100644 --- a/src/qs_wannier90.F +++ b/src/qs_wannier90.F @@ -23,7 +23,6 @@ MODULE qs_wannier90 dbcsr_deallocate_matrix_set USE cp_files, ONLY: close_file,& open_file - USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type USE cp_fm_struct, ONLY: cp_fm_struct_create,& cp_fm_struct_release,& cp_fm_struct_type @@ -53,6 +52,7 @@ MODULE qs_wannier90 dp USE kpoint_methods, ONLY: kpoint_env_initialize,& kpoint_init_cell_index,& + kpoint_initialize_mo_set,& kpoint_initialize_mos,& rskp_transform USE kpoint_types, ONLY: get_kpoint_info,& @@ -70,9 +70,7 @@ MODULE qs_wannier90 qs_env_release,& qs_environment_type USE qs_gamma2kp, ONLY: create_kp_from_gamma - USE qs_matrix_pools, ONLY: mpools_get USE qs_mo_types, ONLY: get_mo_set,& - init_mo_set,& mo_set_p_type USE qs_moments, ONLY: build_berry_kpoint_matrix USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type @@ -160,8 +158,8 @@ SUBROUTINE wannier90_files(qs_env, input, iw) CHARACTER(len=20), ALLOCATABLE, DIMENSION(:) :: atom_symbols CHARACTER(LEN=256) :: datx CHARACTER(len=default_string_length) :: filename, seed_name - INTEGER :: i, i_rep, ib, ib1, ib2, ibs, ic, ik, ik2, ikk, ikpgr, ispin, iunit, ix, iy, iz, & - k, n_rep, nadd, nao, nbs, nexcl, nkp, nmo, nntot, nspins, num_atoms, num_bands, & + INTEGER :: i, i_rep, ib, ib1, ib2, ibs, ik, ik2, ikk, ikpgr, ispin, iunit, ix, iy, iz, k, & + n_rep, nadd, nao, nbs, nexcl, nkp, nmo, nntot, nspins, num_atoms, num_bands, & num_bands_tot, num_kpts, num_wann INTEGER, ALLOCATABLE, DIMENSION(:) :: exclude_bands INTEGER, ALLOCATABLE, DIMENSION(:, :) :: nblist, nnlist @@ -183,10 +181,9 @@ SUBROUTINE wannier90_files(qs_env, input, iw) TYPE(cell_type), POINTER :: cell TYPE(cp_blacs_env_type), POINTER :: blacs_env TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: fmk1, fmk2 - TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools TYPE(cp_fm_struct_type), POINTER :: matrix_struct_mmn, matrix_struct_work TYPE(cp_fm_type), POINTER :: fm_tmp, fmdummy, fmi, fmr, mmn_imag, & - mmn_real, mo_coeff + mmn_real TYPE(cp_para_env_type), POINTER :: para_env TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix @@ -194,7 +191,6 @@ SUBROUTINE wannier90_files(qs_env, input, iw) TYPE(kpoint_env_type), POINTER :: kp TYPE(kpoint_type), POINTER :: kpoint TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos - TYPE(mo_set_p_type), DIMENSION(:, :), POINTER :: moskp TYPE(neighbor_list_set_p_type), DIMENSION(:), & POINTER :: sab_nl TYPE(particle_type), DIMENSION(:), POINTER :: particle_set @@ -356,21 +352,7 @@ SUBROUTINE wannier90_files(qs_env, input, iw) CALL cp_blacs_env_retain(blacs_env) CALL kpoint_env_initialize(kpoint) CALL kpoint_initialize_mos(kpoint, mos, nadd) - DO ik = 1, SIZE(kpoint%kp_env) - CALL mpools_get(kpoint%mpools, ao_mo_fm_pools=ao_mo_fm_pools) - moskp => kpoint%kp_env(ik)%kpoint_env%mos - ikk = kpoint%kp_range(1) + ik - 1 - CPASSERT(ASSOCIATED(moskp)) - DO ispin = 1, SIZE(moskp, 2) - DO ic = 1, SIZE(moskp, 1) - CALL get_mo_set(moskp(ic, ispin)%mo_set, mo_coeff=mo_coeff) - IF (.NOT. ASSOCIATED(mo_coeff)) THEN - CALL init_mo_set(moskp(ic, ispin)%mo_set, & - fm_pool=ao_mo_fm_pools(ispin)%pool, name="kpoints") - END IF - END DO - END DO - END DO + CALL kpoint_initialize_mo_set(kpoint) ! CALL get_qs_env(qs_env=qs_env_kp, sab_orb=sab_nl, dft_control=dft_control) CALL kpoint_init_cell_index(kpoint, sab_nl, para_env, dft_control) diff --git a/src/rpa_gw.F b/src/rpa_gw.F index 4bbcfbfbeb..794b36f059 100644 --- a/src/rpa_gw.F +++ b/src/rpa_gw.F @@ -1381,14 +1381,14 @@ SUBROUTINE compute_QP_energies(vec_Sigma_c_gw, count_ev_sc_GW, gw_corr_lev_occ, CALL print_and_update_for_ev_sc( & vec_gw_energ(:, 1), & z_value(:, 1), m_value(:, 1), mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp), Eigenval(:, 1), & - Eigenval_last(:, 1), Eigenval_scf(:, 1), gw_corr_lev_occ(1), gw_corr_lev_tot, & + Eigenval_last(:, 1), Eigenval_scf(:, 1), gw_corr_lev_occ(1), gw_corr_lev_virt(1), gw_corr_lev_tot, & crossing_search, homo(1), unit_nr, count_ev_sc_GW_print, count_sc_GW0_print, & ikp, nkp_self_energy, kpoints, 1) CALL print_and_update_for_ev_sc( & vec_gw_energ(:, 2), & z_value(:, 2), m_value(:, 2), mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 2, ikp), Eigenval(:, 2), & - Eigenval_last(:, 2), Eigenval_scf(:, 2), gw_corr_lev_occ(2), gw_corr_lev_tot, & + Eigenval_last(:, 2), Eigenval_scf(:, 2), gw_corr_lev_occ(2), gw_corr_lev_virt(2), gw_corr_lev_tot, & crossing_search, homo(2), unit_nr, count_ev_sc_GW_print, count_sc_GW0_print, & ikp, nkp_self_energy, kpoints, 2) @@ -1409,7 +1409,7 @@ SUBROUTINE compute_QP_energies(vec_Sigma_c_gw, count_ev_sc_GW, gw_corr_lev_occ, CALL print_and_update_for_ev_sc( & vec_gw_energ(:, 1), & z_value(:, 1), m_value(:, 1), mp2_env%ri_g0w0%vec_Sigma_x_minus_vxc_gw(:, 1, ikp), Eigenval(:, 1), & - Eigenval_last(:, 1), Eigenval_scf(:, 1), gw_corr_lev_occ(1), gw_corr_lev_tot, & + Eigenval_last(:, 1), Eigenval_scf(:, 1), gw_corr_lev_occ(1), gw_corr_lev_virt(1), gw_corr_lev_tot, & crossing_search, homo(1), unit_nr, count_ev_sc_GW_print, count_sc_GW0_print, & ikp, nkp_self_energy, kpoints, 0) @@ -3765,6 +3765,7 @@ END SUBROUTINE get_sigma_c_newton_pade !> \param Eigenval_last ... !> \param Eigenval_scf ... !> \param gw_corr_lev_occ ... +!> \param gw_corr_lev_virt ... !> \param gw_corr_lev_tot ... !> \param crossing_search ... !> \param homo ... @@ -3778,15 +3779,16 @@ END SUBROUTINE get_sigma_c_newton_pade ! ************************************************************************************************** SUBROUTINE print_and_update_for_ev_sc(vec_gw_energ, & z_value, m_value, vec_Sigma_x_minus_vxc_gw, Eigenval, & - Eigenval_last, Eigenval_scf, gw_corr_lev_occ, gw_corr_lev_tot, & + Eigenval_last, Eigenval_scf, & + gw_corr_lev_occ, gw_corr_lev_virt, gw_corr_lev_tot, & crossing_search, homo, unit_nr, count_ev_sc_GW, count_sc_GW0, & ikp, nkp_self_energy, kpoints, ispin) REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vec_gw_energ, z_value, m_value REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: vec_Sigma_x_minus_vxc_gw, Eigenval, & Eigenval_last, Eigenval_scf - INTEGER, INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_tot, crossing_search, homo, unit_nr, & - count_ev_sc_GW, count_sc_GW0, ikp, nkp_self_energy + INTEGER, INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, gw_corr_lev_tot, crossing_search, & + homo, unit_nr, count_ev_sc_GW, count_sc_GW0, ikp, nkp_self_energy TYPE(kpoint_type), INTENT(IN), POINTER :: kpoints INTEGER, INTENT(IN) :: ispin @@ -3796,7 +3798,8 @@ SUBROUTINE print_and_update_for_ev_sc(vec_gw_energ, & INTEGER :: handle, n_level_gw, n_level_gw_ref LOGICAL :: do_alpha, do_beta, do_closed_shell, & do_kpoints, is_energy_okay - REAL(KIND=dp) :: new_energy + REAL(KIND=dp) :: E_GAP_GW, E_HOMO_GW, E_LUMO_GW, & + new_energy CALL timeset(routineN, handle) @@ -3937,15 +3940,19 @@ SUBROUTINE print_and_update_for_ev_sc(vec_gw_energ, & END DO IF (unit_nr > 0) THEN + E_HOMO_GW = MAXVAL(Eigenval(homo - gw_corr_lev_occ + 1:homo)) + E_LUMO_GW = MINVAL(Eigenval(homo + 1:homo + gw_corr_lev_virt)) + E_GAP_GW = E_LUMO_GW - E_HOMO_GW + IF (do_closed_shell) THEN WRITE (unit_nr, '(T3,A)') ' ' - WRITE (unit_nr, '(T3,A,F57.2)') 'GW HOMO-LUMO gap (eV)', (Eigenval(homo + 1) - Eigenval(homo))*evolt + WRITE (unit_nr, '(T3,A,F57.2)') 'GW HOMO-LUMO gap (eV)', E_GAP_GW*evolt ELSE IF (do_alpha) THEN WRITE (unit_nr, '(T3,A)') ' ' - WRITE (unit_nr, '(T3,A,F51.2)') 'Alpha GW HOMO-LUMO gap (eV)', (Eigenval(homo + 1) - Eigenval(homo))*evolt + WRITE (unit_nr, '(T3,A,F51.2)') 'Alpha GW HOMO-LUMO gap (eV)', E_GAP_GW*evolt ELSE IF (do_beta) THEN WRITE (unit_nr, '(T3,A)') ' ' - WRITE (unit_nr, '(T3,A,F52.2)') 'Beta GW HOMO-LUMO gap (eV)', (Eigenval(homo + 1) - Eigenval(homo))*evolt + WRITE (unit_nr, '(T3,A,F52.2)') 'Beta GW HOMO-LUMO gap (eV)', E_GAP_GW*evolt END IF END IF diff --git a/src/rpa_gw_im_time_util.F b/src/rpa_gw_im_time_util.F index e17bbae861..f9bd787484 100644 --- a/src/rpa_gw_im_time_util.F +++ b/src/rpa_gw_im_time_util.F @@ -41,6 +41,8 @@ MODULE rpa_gw_im_time_util hfx_compression_type USE kinds, ONLY: dp,& int_8 + USE mathconstants, ONLY: pi,& + twopi USE message_passing, ONLY: mp_alltoall,& mp_dims_create,& mp_sum @@ -67,7 +69,8 @@ MODULE rpa_gw_im_time_util CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_gw_im_time_util' - PUBLIC :: get_tensor_3c_overl_int_gw + PUBLIC :: get_tensor_3c_overl_int_gw, compute_weight_re_im, get_atom_index_from_basis_function_index, & + compute_shortest_lattice_vector CONTAINS ! ************************************************************************************************** @@ -417,7 +420,7 @@ SUBROUTINE get_tensor_3c_overl_int_gw(t_3c_overl_int, t_3c_O_compressed, t_3c_O_ IF (unit_nr > 0) THEN WRITE (UNIT=unit_nr, FMT="((T3,A,T66,F11.2,A4))") & - "MEMORY_INFO| Memory for MO-contracted integral tensor (compressed):", memory_3c, ' MiB' + "MEMORY_INFO| Memory of MO-contracted tensor (compressed):", memory_3c, ' MiB' WRITE (UNIT=unit_nr, FMT="((T3,A,T60,F21.2))") & "MEMORY_INFO| Compression factor: ", compression_factor @@ -733,4 +736,174 @@ SUBROUTINE reflect_mat_row(mat_reflected, mat_orig, para_env, qs_env, unit_nr, d END SUBROUTINE +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param atom_from_basis_index ... +!> \param basis_size ... +!> \param basis_type ... +! ************************************************************************************************** + SUBROUTINE get_atom_index_from_basis_function_index(qs_env, atom_from_basis_index, basis_size, basis_type) + TYPE(qs_environment_type), POINTER :: qs_env + INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_basis_index + INTEGER :: basis_size + CHARACTER(LEN=*) :: basis_type + + CHARACTER(LEN=*), PARAMETER :: routineN = 'get_atom_index_from_basis_function_index' + + INTEGER :: iatom, LLL, natom, nkind + INTEGER, DIMENSION(:), POINTER :: row_blk_end, row_blk_start + TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + + NULLIFY (qs_kind_set, particle_set) + CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, natom=natom, nkind=nkind, & + particle_set=particle_set) + + ALLOCATE (row_blk_start(natom)) + ALLOCATE (row_blk_end(natom)) + ALLOCATE (basis_set(nkind)) + CALL basis_set_list_setup(basis_set, basis_type, qs_kind_set) + CALL get_particle_set(particle_set, qs_kind_set, first_sgf=row_blk_start, last_sgf=row_blk_end, & + basis=basis_set) + DO LLL = 1, basis_size + DO iatom = 1, natom + IF (LLL >= row_blk_start(iatom) .AND. LLL <= row_blk_end(iatom)) THEN + atom_from_basis_index(LLL) = iatom + END IF + END DO + END DO + + DEALLOCATE (basis_set) + DEALLOCATE (row_blk_start) + DEALLOCATE (row_blk_end) + + END SUBROUTINE get_atom_index_from_basis_function_index + +! ************************************************************************************************** +!> \brief ... +!> \param weight_re ... +!> \param weight_im ... +!> \param num_cells ... +!> \param iatom ... +!> \param jatom ... +!> \param xkp ... +!> \param wkp_W ... +!> \param cell ... +!> \param index_to_cell ... +!> \param hmat ... +!> \param particle_set ... +!> \param abs_cutoffs_chi_W ... +! ************************************************************************************************** + SUBROUTINE compute_weight_re_im(weight_re, weight_im, & + num_cells, iatom, jatom, xkp, wkp_W, & + cell, index_to_cell, hmat, particle_set, abs_cutoffs_chi_W) + + REAL(KIND=dp) :: weight_re, weight_im + INTEGER :: num_cells, iatom, jatom + REAL(KIND=dp), DIMENSION(3) :: xkp + REAL(KIND=dp) :: wkp_W + TYPE(cell_type), POINTER :: cell + INTEGER, DIMENSION(:, :), POINTER :: index_to_cell + REAL(KIND=dp), DIMENSION(3, 3) :: hmat + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + REAL(KIND=dp), DIMENSION(2) :: abs_cutoffs_chi_W + + CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_weight_re_im' + + INTEGER :: handle, icell, xcell, ycell, zcell + REAL(KIND=dp) :: abs_rab_cell, arg, coskl, & + dampening_function, first_cutoff, & + second_cutoff, sinkl + REAL(KIND=dp), DIMENSION(3) :: cell_vector, rab_cell_i + + CALL timeset(routineN, handle) + + first_cutoff = abs_cutoffs_chi_W(1) + second_cutoff = abs_cutoffs_chi_W(2) + + weight_re = 0.0_dp + weight_im = 0.0_dp + + DO icell = 1, num_cells + + xcell = index_to_cell(1, icell) + ycell = index_to_cell(2, icell) + zcell = index_to_cell(3, icell) + + arg = REAL(xcell, dp)*xkp(1) + REAL(ycell, dp)*xkp(2) + REAL(zcell, dp)*xkp(3) + + coskl = wkp_W*COS(twopi*arg) + sinkl = wkp_W*SIN(twopi*arg) + + cell_vector(1:3) = MATMUL(hmat, REAL(index_to_cell(1:3, icell), dp)) + + rab_cell_i(1:3) = pbc(particle_set(iatom)%r(1:3), cell) - & + (pbc(particle_set(jatom)%r(1:3), cell) + cell_vector(1:3)) + + abs_rab_cell = SQRT(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2) + + ! smooth function that is only one if atoms are very close and becomes smoothly zero + ! for larger distances + IF (abs_rab_cell < first_cutoff) THEN + weight_re = 1.0_dp*coskl + weight_im = 1.0_dp*sinkl + ELSE IF (abs_rab_cell < second_cutoff) THEN + dampening_function = 0.5_dp*(COS(pi*(abs_rab_cell - first_cutoff)/ & + (second_cutoff - first_cutoff)) + 1.0_dp) + weight_re = dampening_function*coskl + weight_im = dampening_function*sinkl + END IF + + END DO + + CALL timestop(handle) + + END SUBROUTINE compute_weight_re_im + +! ************************************************************************************************** +!> \brief ... +!> \param shortest_lattice_vector ... +!> \param hmat ... +!> \param index_to_cell ... +!> \param num_cells ... +! ************************************************************************************************** + SUBROUTINE compute_shortest_lattice_vector(shortest_lattice_vector, hmat, index_to_cell, num_cells) + + REAL(KIND=dp) :: shortest_lattice_vector + REAL(KIND=dp), DIMENSION(3, 3) :: hmat + INTEGER, DIMENSION(:, :), POINTER :: index_to_cell + INTEGER :: num_cells + + CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_shortest_lattice_vector' + + INTEGER :: handle, icell, xcell, ycell, zcell + REAL(KIND=dp) :: abs_cell_vector + REAL(KIND=dp), DIMENSION(3) :: cell_vector + + CALL timeset(routineN, handle) + + shortest_lattice_vector = 1.0E9_dp + + DO icell = 1, num_cells + + xcell = index_to_cell(1, icell) + ycell = index_to_cell(2, icell) + zcell = index_to_cell(3, icell) + + IF (xcell == 0 .AND. ycell == 0 .AND. zcell == 0) CYCLE + + cell_vector(1:3) = MATMUL(hmat, REAL(index_to_cell(1:3, icell), dp)) + + abs_cell_vector = SQRT(cell_vector(1)**2 + cell_vector(2)**2 + cell_vector(3)**2) + + IF (abs_cell_vector < shortest_lattice_vector) shortest_lattice_vector = abs_cell_vector + + END DO + + CALL timestop(handle) + + END SUBROUTINE compute_shortest_lattice_vector + END MODULE rpa_gw_im_time_util diff --git a/src/rpa_gw_kpoints.F b/src/rpa_gw_kpoints.F index 124d257038..c506a51e0b 100644 --- a/src/rpa_gw_kpoints.F +++ b/src/rpa_gw_kpoints.F @@ -11,10 +11,8 @@ !> 04.2019 created [Jan Wilhelm] ! ************************************************************************************************** MODULE rpa_gw_kpoints - USE basis_set_types, ONLY: gto_basis_set_p_type USE cell_types, ONLY: cell_type,& - get_cell,& - pbc + get_cell USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_invert,& cp_cfm_gemm,& cp_cfm_scale_and_add,& @@ -39,20 +37,21 @@ MODULE rpa_gw_kpoints dbcsr_create,& dbcsr_p_type,& dbcsr_release + USE input_constants, ONLY: kp_weights_W_auto,& + kp_weights_W_tailored,& + kp_weights_W_uniform USE kinds, ONLY: dp USE kpoint_types, ONLY: get_kpoint_info,& kpoint_type USE mathconstants, ONLY: gaussi,& - twopi,& z_one,& z_zero USE mathlib, ONLY: invmat - USE particle_methods, ONLY: get_particle_set USE particle_types, ONLY: particle_type USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type - USE qs_integral_utils, ONLY: basis_set_list_setup - USE qs_kind_types, ONLY: qs_kind_type + USE rpa_gw_im_time_util, ONLY: compute_weight_re_im,& + get_atom_index_from_basis_function_index #include "./base/base_uses.f90" IMPLICIT NONE @@ -110,17 +109,15 @@ SUBROUTINE compute_Wc_real_space_tau_GW(fm_mat_W_tau, cfm_mat_Q, fm_mat_L_re, fm CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Wc_real_space_tau_GW' - INTEGER :: handle, handle2, i_global, iatom, iatom_old, icell, iiB, iquad, irow, j_global, & - jatom, jatom_old, jcol, jjB, jkp, LLL, natom, ncol_local, nkind, nkp, nrow_local, & - num_cells, xcell, ycell, zcell + INTEGER :: handle, handle2, i_global, iatom, iatom_old, iiB, iquad, irow, j_global, jatom, & + jatom_old, jcol, jjB, jkp, ncol_local, nkp, nrow_local, num_cells INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_RI_index - INTEGER, DIMENSION(:), POINTER :: col_indices, row_blk_end, row_blk_start, & - row_indices + INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices INTEGER, DIMENSION(:, :), POINTER :: index_to_cell LOGICAL :: do_V_and_not_W - REAL(KIND=dp) :: abs_rab_cell, arg, contribution, coskl, cutoff_exp, d_0, omega, sinkl, & - sum_exp, sum_exp_k_im, sum_exp_k_re, tau, weight, weight_im, weight_re - REAL(KIND=dp), DIMENSION(3) :: cell_vector, rab_cell_i + REAL(KIND=dp) :: contribution, omega, tau, weight, & + weight_im, weight_re + REAL(KIND=dp), DIMENSION(2) :: abs_cutoffs_chi_W REAL(KIND=dp), DIMENSION(3, 3) :: hmat REAL(KIND=dp), DIMENSION(:), POINTER :: wkp REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp @@ -128,9 +125,7 @@ SUBROUTINE compute_Wc_real_space_tau_GW(fm_mat_W_tau, cfm_mat_Q, fm_mat_L_re, fm TYPE(cp_cfm_type), POINTER :: cfm_mat_L, cfm_mat_work, cfm_mat_work_2 TYPE(cp_fm_type), POINTER :: fm_dummy, fm_mat_work_global, & fm_mat_work_local - TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_RI_tmp TYPE(particle_type), DIMENSION(:), POINTER :: particle_set - TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set CALL timeset(routineN, handle) @@ -216,29 +211,16 @@ SUBROUTINE compute_Wc_real_space_tau_GW(fm_mat_W_tau, cfm_mat_Q, fm_mat_L_re, fm CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp, nkp=nkp) index_to_cell => kpoints%index_to_cell num_cells = SIZE(index_to_cell, 2) - d_0 = qs_env%mp2_env%ri_rpa_im_time%cutoff - cutoff_exp = 10000.0_dp + + abs_cutoffs_chi_W = qs_env%mp2_env%ri_rpa_im_time%abs_cutoffs_chi_W CALL cp_cfm_set_all(cfm_mat_work, z_zero) - NULLIFY (qs_kind_set, cell, particle_set) - CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, cell=cell, natom=natom, nkind=nkind, & - particle_set=particle_set) - - ALLOCATE (row_blk_start(natom)) - ALLOCATE (row_blk_end(natom)) - ALLOCATE (basis_set_RI_tmp(nkind)) - CALL basis_set_list_setup(basis_set_RI_tmp, "RI_AUX", qs_kind_set) - CALL get_particle_set(particle_set, qs_kind_set, first_sgf=row_blk_start, last_sgf=row_blk_end, & - basis=basis_set_RI_tmp) - DEALLOCATE (basis_set_RI_tmp) ALLOCATE (atom_from_RI_index(dimen_RI)) - DO LLL = 1, dimen_RI - DO iatom = 1, natom - IF (LLL >= row_blk_start(iatom) .AND. LLL <= row_blk_end(iatom)) THEN - atom_from_RI_index(LLL) = iatom - END IF - END DO - END DO + + CALL get_atom_index_from_basis_function_index(qs_env, atom_from_RI_index, dimen_RI, "RI_AUX") + + NULLIFY (cell, particle_set) + CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set) CALL get_cell(cell=cell, h=hmat) iatom_old = 0 jatom_old = 0 @@ -257,38 +239,9 @@ SUBROUTINE compute_Wc_real_space_tau_GW(fm_mat_W_tau, cfm_mat_Q, fm_mat_L_re, fm IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old) THEN - sum_exp = 0.0_dp - sum_exp_k_re = 0.0_dp - sum_exp_k_im = 0.0_dp - - DO icell = 1, num_cells - - xcell = index_to_cell(1, icell) - ycell = index_to_cell(2, icell) - zcell = index_to_cell(3, icell) - - arg = REAL(xcell, dp)*xkp(1, ikp) + REAL(ycell, dp)*xkp(2, ikp) + REAL(zcell, dp)*xkp(3, ikp) - - coskl = wkp_W(ikp)*COS(twopi*arg) - sinkl = wkp_W(ikp)*SIN(twopi*arg) - - cell_vector(1:3) = MATMUL(hmat, REAL(index_to_cell(1:3, icell), dp)) - - rab_cell_i(1:3) = pbc(particle_set(iatom)%r(1:3), cell) - & - (pbc(particle_set(jatom)%r(1:3), cell) + cell_vector(1:3)) - - abs_rab_cell = SQRT(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2) - - IF (abs_rab_cell/d_0 < cutoff_exp) THEN - sum_exp = sum_exp + EXP(-abs_rab_cell/d_0) - sum_exp_k_re = sum_exp_k_re + EXP(-abs_rab_cell/d_0)*coskl - sum_exp_k_im = sum_exp_k_im + EXP(-abs_rab_cell/d_0)*sinkl - END IF - - END DO - - weight_re = sum_exp_k_re/sum_exp - weight_im = sum_exp_k_im/sum_exp + CALL compute_weight_re_im(weight_re, weight_im, & + num_cells, iatom, jatom, xkp(1:3, ikp), wkp_W(ikp), & + cell, index_to_cell, hmat, particle_set, abs_cutoffs_chi_W) iatom_old = iatom jatom_old = jatom @@ -384,14 +337,12 @@ SUBROUTINE compute_Wc_real_space_tau_GW(fm_mat_W_tau, cfm_mat_Q, fm_mat_L_re, fm CALL cp_fm_release(fm_mat_work_global) CALL cp_fm_release(fm_mat_work_local) DEALLOCATE (atom_from_RI_index) - DEALLOCATE (row_blk_start) - DEALLOCATE (row_blk_end) CALL timestop(handle2) CALL timestop(handle) - END SUBROUTINE + END SUBROUTINE compute_Wc_real_space_tau_GW ! ************************************************************************************************** !> \brief ... @@ -526,33 +477,36 @@ SUBROUTINE compute_Wc_kp_tau_GW(cfm_mat_W_kp_tau, cfm_mat_Q, fm_mat_L_re, fm_mat ! ************************************************************************************************** !> \brief ... +!> \param qs_env ... !> \param wkp_W ... +!> \param wkp_V ... !> \param kpoints ... -!> \param h_mat ... !> \param h_inv ... -!> \param exp_kpoints ... !> \param periodic ... ! ************************************************************************************************** - SUBROUTINE compute_wkp_W(wkp_W, kpoints, h_mat, h_inv, exp_kpoints, periodic) + SUBROUTINE compute_wkp_W(qs_env, wkp_W, wkp_V, kpoints, h_inv, periodic) + + TYPE(qs_environment_type), POINTER :: qs_env REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & - INTENT(OUT) :: wkp_W - TYPE(kpoint_type), INTENT(IN), POINTER :: kpoints - REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: h_mat, h_inv - REAL(KIND=dp), INTENT(INOUT) :: exp_kpoints - INTEGER, DIMENSION(3), INTENT(IN) :: periodic + INTENT(OUT) :: wkp_W, wkp_V + TYPE(kpoint_type), POINTER :: kpoints + REAL(KIND=dp), DIMENSION(3, 3) :: h_inv + INTEGER, DIMENSION(3) :: periodic CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_wkp_W' - INTEGER :: handle, i_dim, i_x, ikp, info, j_y, k_z, & - n_x, n_y, n_z, nkp, nsuperfine, & - num_lin_eqs - REAL(KIND=dp) :: a_vec_dot_k_vec, integral, k_sq, weight - REAL(KIND=dp), DIMENSION(3) :: a_vec, k_vec, x_vec - REAL(KIND=dp), DIMENSION(:), POINTER :: right_side, wkp + INTEGER :: handle, i_x, ikp, info, j_y, k_z, & + kpoint_weights_W_method, n_x, n_y, & + n_z, nkp, nsuperfine, num_lin_eqs + REAL(KIND=dp) :: exp_kpoints, integral, k_sq, weight + REAL(KIND=dp), DIMENSION(3) :: k_vec, x_vec + REAL(KIND=dp), DIMENSION(:), POINTER :: right_side, wkp, wkp_tmp REAL(KIND=dp), DIMENSION(:, :), POINTER :: matrix_lin_eqs, xkp CALL timeset(routineN, handle) + kpoint_weights_W_method = qs_env%mp2_env%ri_rpa_im_time%kpoint_weights_W_method + CALL get_kpoint_info(kpoints, xkp=xkp, wkp=wkp, nkp=nkp) ! we determine the kpoint weights of the Monkhors Pack mesh new @@ -563,158 +517,120 @@ SUBROUTINE compute_wkp_W(wkp_W, kpoints, h_mat, h_inv, exp_kpoints, periodic) ! 1) 1/k^2, 1/k and const are integrated exactly ! 2) the kpoint weights of kpoints with identical absolute value are ! the same, of e.g. (1/8,3/8,3/8) same weight as for (3/8,1/8,3/8) - ! for 1d and 2d materials: we use normal Monkhorst-Pack mesh, checked + ! for 1d and 2d materials: we use ordinary Monkhorst-Pack weights, checked ! by SUM(periodic) == 3 + ALLOCATE (wkp_V(nkp), wkp_W(nkp)) - IF (exp_kpoints < 2.0_dp .AND. SUM(periodic) == 3) THEN - - ! first, compute the integral of f(k)=1/k^2 and 1/k on super fine grid - nsuperfine = 500 - integral = 0.0_dp - IF (exp_kpoints > 0.0_dp) exp_kpoints = -2.0_dp - - ! actually, there is the factor *det_3x3(h_inv) missing to account for the - ! integration volume but for wkp det_3x3(h_inv) is needed - weight = 2.0_dp/(REAL(nsuperfine, dp))**3 - DO i_x = 1, nsuperfine - DO j_y = 1, nsuperfine - DO k_z = 1, nsuperfine/2 - - x_vec = (/REAL(i_x - nsuperfine/2, dp) - 0.5_dp, & - REAL(j_y - nsuperfine/2, dp) - 0.5_dp, & - REAL(k_z - nsuperfine/2, dp) - 0.5_dp/)/ & - REAL(nsuperfine, dp) - k_vec = MATMUL(h_inv(1:3, 1:3), x_vec) - k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2 - integral = integral + weight*k_sq**(exp_kpoints*0.5_dp) - END DO - END DO - END DO - - num_lin_eqs = nkp + 2 - - ALLOCATE (matrix_lin_eqs(num_lin_eqs, num_lin_eqs)) - matrix_lin_eqs(:, :) = 0.0_dp + IF (kpoint_weights_W_method == kp_weights_W_uniform) THEN - DO ikp = 1, nkp - - k_vec = MATMUL(h_inv(1:3, 1:3), xkp(1:3, ikp)) - k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2 - - matrix_lin_eqs(ikp, ikp) = 2.0_dp - matrix_lin_eqs(ikp, nkp + 1) = 1.0_dp - matrix_lin_eqs(ikp, nkp + 2) = 1.0_dp*k_sq**(exp_kpoints*0.5_dp) - - matrix_lin_eqs(nkp + 1, ikp) = 1.0_dp - matrix_lin_eqs(nkp + 2, ikp) = 1.0_dp*k_sq**(exp_kpoints*0.5_dp) - - END DO - - CALL invmat(matrix_lin_eqs, info) - ! check whether inversion was successful - CPASSERT(info == 0) - - ALLOCATE (wkp_W(num_lin_eqs)) + wkp_V(:) = wkp(:) + wkp_W(:) = wkp(:) - ALLOCATE (right_side(num_lin_eqs)) - right_side = 0.0_dp - right_side(nkp + 1) = 1.0_dp - right_side(nkp + 2) = integral + ELSE IF (kpoint_weights_W_method == kp_weights_W_tailored .OR. & + kpoint_weights_W_method == kp_weights_W_auto) THEN - wkp_W(1:num_lin_eqs) = MATMUL(matrix_lin_eqs, right_side) + IF (SUM(periodic) == 1 .OR. SUM(periodic) == 2) THEN - DEALLOCATE (matrix_lin_eqs, right_side) + wkp_V(:) = wkp(:) + wkp_W(:) = wkp(:) - ELSE IF (exp_kpoints < 2.0_dp .AND. SUM(periodic) == 1) THEN + ELSE IF (SUM(periodic) == 3) THEN - ! first, compute the integral of f(k)=1/k^2 and 1/k on super fine grid - nsuperfine = 5000 - integral = 0.0_dp + IF (kpoint_weights_W_method == kp_weights_W_tailored) & + exp_kpoints = qs_env%mp2_env%ri_rpa_im_time%exp_tailored_weights + IF (kpoint_weights_W_method == kp_weights_W_auto) exp_kpoints = -2.0_dp - ! actually, there is the factor *det_3x3(h_inv) missing to account for the - ! integration volume but for wkp det_3x3(h_inv) is needed - weight = 1.0_dp/REAL(nsuperfine, dp) - IF (periodic(1) == 1) THEN - n_x = nsuperfine - ELSE - n_x = 1 - END IF - IF (periodic(2) == 1) THEN - n_y = nsuperfine - ELSE - n_y = 1 - END IF - IF (periodic(3) == 1) THEN - n_z = nsuperfine - ELSE - n_z = 1 - END IF + ! first, compute the integral of f(k)=1/k^2 and 1/k on super fine grid + nsuperfine = 500 + integral = 0.0_dp - a_vec = MATMUL(h_mat(1:3, 1:3), & - (/REAL(periodic(1), dp), REAL(periodic(2), dp), REAL(periodic(3), dp)/)) + IF (periodic(1) == 1) THEN + n_x = nsuperfine + ELSE + n_x = 1 + END IF + IF (periodic(2) == 1) THEN + n_y = nsuperfine + ELSE + n_y = 1 + END IF + IF (periodic(3) == 1) THEN + n_z = nsuperfine + ELSE + n_z = 1 + END IF - DO i_x = 1, n_x - DO j_y = 1, n_y - DO k_z = 1, n_z + ! actually, there is the factor *det_3x3(h_inv) missing to account for the + ! integration volume but for wkp det_3x3(h_inv) is needed + weight = 1.0_dp/(REAL(n_x, dp)*REAL(n_y, dp)*REAL(n_z, dp)) + DO i_x = 1, n_x + DO j_y = 1, n_y + DO k_z = 1, n_z + + IF (periodic(1) == 1) THEN + x_vec(1) = (REAL(i_x - nsuperfine/2, dp) - 0.5_dp)/REAL(nsuperfine, dp) + ELSE + x_vec(1) = 0.0_dp + END IF + IF (periodic(2) == 1) THEN + x_vec(2) = (REAL(j_y - nsuperfine/2, dp) - 0.5_dp)/REAL(nsuperfine, dp) + ELSE + x_vec(2) = 0.0_dp + END IF + IF (periodic(3) == 1) THEN + x_vec(3) = (REAL(k_z - nsuperfine/2, dp) - 0.5_dp)/REAL(nsuperfine, dp) + ELSE + x_vec(3) = 0.0_dp + END IF - x_vec = (/REAL(i_x - nsuperfine/2, dp) - 0.5_dp, & - REAL(j_y - nsuperfine/2, dp) - 0.5_dp, & - REAL(k_z - nsuperfine/2, dp) - 0.5_dp/)/ & - REAL(nsuperfine, dp) + k_vec = MATMUL(h_inv(1:3, 1:3), x_vec) + k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2 + integral = integral + weight*k_sq**(exp_kpoints*0.5_dp) - DO i_dim = 1, 3 - IF (periodic(i_dim) == 0) THEN - x_vec(i_dim) = 0.0_dp - END IF END DO - - k_vec = MATMUL(h_inv(1:3, 1:3), x_vec) - a_vec_dot_k_vec = a_vec(1)*k_vec(1) + a_vec(2)*k_vec(2) + a_vec(3)*k_vec(3) - integral = integral + weight*LOG(2.0_dp - 2.0_dp*COS(a_vec_dot_k_vec)) END DO END DO - END DO - num_lin_eqs = nkp + 2 + num_lin_eqs = nkp + 2 - ALLOCATE (matrix_lin_eqs(num_lin_eqs, num_lin_eqs)) - matrix_lin_eqs(:, :) = 0.0_dp + ALLOCATE (matrix_lin_eqs(num_lin_eqs, num_lin_eqs)) + matrix_lin_eqs(:, :) = 0.0_dp - DO ikp = 1, nkp + DO ikp = 1, nkp - k_vec = MATMUL(h_inv(1:3, 1:3), xkp(1:3, ikp)) - k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2 + k_vec = MATMUL(h_inv(1:3, 1:3), xkp(1:3, ikp)) + k_sq = k_vec(1)**2 + k_vec(2)**2 + k_vec(3)**2 - matrix_lin_eqs(ikp, ikp) = 2.0_dp - matrix_lin_eqs(ikp, nkp + 1) = 1.0_dp + matrix_lin_eqs(ikp, ikp) = 2.0_dp + matrix_lin_eqs(ikp, nkp + 1) = 1.0_dp + matrix_lin_eqs(nkp + 1, ikp) = 1.0_dp - a_vec_dot_k_vec = a_vec(1)*k_vec(1) + a_vec(2)*k_vec(2) + a_vec(3)*k_vec(3) - matrix_lin_eqs(ikp, nkp + 2) = LOG(2.0_dp - 2.0_dp*COS(a_vec_dot_k_vec)) + matrix_lin_eqs(ikp, nkp + 2) = k_sq**(exp_kpoints*0.5_dp) + matrix_lin_eqs(nkp + 2, ikp) = k_sq**(exp_kpoints*0.5_dp) - matrix_lin_eqs(nkp + 1, ikp) = 1.0_dp - matrix_lin_eqs(nkp + 2, ikp) = LOG(2.0_dp - 2.0_dp*COS(a_vec_dot_k_vec)) + END DO - END DO + CALL invmat(matrix_lin_eqs, info) + ! check whether inversion was successful + CPASSERT(info == 0) - CALL invmat(matrix_lin_eqs, info) - ! check whether inversion was successful - CPASSERT(info == 0) + ALLOCATE (right_side(num_lin_eqs)) + right_side = 0.0_dp + right_side(nkp + 1) = 1.0_dp + ! divide integral by two because CP2K k-mesh already considers symmetry k <-> -k + right_side(nkp + 2) = integral - ALLOCATE (wkp_W(num_lin_eqs)) + ALLOCATE (wkp_tmp(num_lin_eqs)) - ALLOCATE (right_side(num_lin_eqs)) - right_side = 0.0_dp - right_side(nkp + 1) = 1.0_dp - right_side(nkp + 2) = integral + wkp_tmp(1:num_lin_eqs) = MATMUL(matrix_lin_eqs, right_side) - wkp_W(1:num_lin_eqs) = MATMUL(matrix_lin_eqs, right_side) + wkp_V(1:nkp) = wkp_tmp(1:nkp) - DEALLOCATE (matrix_lin_eqs, right_side) + DEALLOCATE (matrix_lin_eqs, right_side, wkp_tmp) - ELSE + wkp_W(:) = wkp_V(:) - ALLOCATE (wkp_W(nkp)) - wkp_W(:) = wkp(:) + END IF END IF diff --git a/src/rpa_im_time.F b/src/rpa_im_time.F index aad564f183..e4682e8a61 100644 --- a/src/rpa_im_time.F +++ b/src/rpa_im_time.F @@ -6,15 +6,13 @@ !--------------------------------------------------------------------------------------------------! ! ************************************************************************************************** -!> \brief Routines for RPA with imaginary time +!> \brief Routines for low-scaling RPA/GW with imaginary time !> \par History !> 10.2015 created [Jan Wilhelm] ! ************************************************************************************************** MODULE rpa_im_time - USE cell_types, ONLY: cell_type,& get_cell - USE cp_control_types, ONLY: dft_control_type USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,& dbcsr_allocate_matrix_set,& dbcsr_deallocate_matrix_set @@ -24,6 +22,7 @@ MODULE rpa_im_time USE cp_fm_types, ONLY: cp_fm_create,& cp_fm_get_info,& cp_fm_release,& + cp_fm_set_all,& cp_fm_to_fm,& cp_fm_type USE cp_gemm_interface, ONLY: cp_gemm @@ -49,15 +48,18 @@ MODULE rpa_im_time USE mathconstants, ONLY: twopi USE message_passing, ONLY: mp_sync USE mp2_types, ONLY: mp2_type + USE particle_types, ONLY: particle_type USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type USE qs_mo_types, ONLY: get_mo_set,& mo_set_p_type,& mo_set_type - USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type USE qs_tensors, ONLY: decompress_tensor,& get_tensor_occupancy USE qs_tensors_types, ONLY: create_2c_tensor + USE rpa_gw_im_time_util, ONLY: compute_shortest_lattice_vector,& + compute_weight_re_im,& + get_atom_index_from_basis_function_index #include "./base/base_uses.f90" IMPLICIT NONE @@ -108,6 +110,7 @@ MODULE rpa_im_time !> \param mp2_env ... !> \param para_env ... !> \param qs_env ... +!> \param do_kpoints_from_Gamma ... !> \param index_to_cell_3c ... !> \param cell_to_index_3c ... !> \param has_mat_P_blocks ... @@ -129,7 +132,8 @@ SUBROUTINE compute_mat_P_omega(mat_P_omega, fm_scaled_dm_occ_tau, & alpha, eps_filter_im_time, Eigenval, nmo, & num_integ_points, cut_memory, unit_nr, & mp2_env, para_env, & - qs_env, index_to_cell_3c, cell_to_index_3c, & + qs_env, do_kpoints_from_Gamma, & + index_to_cell_3c, cell_to_index_3c, & has_mat_P_blocks, do_ri_sos_laplace_mp2, & dbcsr_time, dbcsr_nflop) TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN) :: mat_P_omega @@ -161,6 +165,7 @@ SUBROUTINE compute_mat_P_omega(mat_P_omega, fm_scaled_dm_occ_tau, & TYPE(mp2_type), POINTER :: mp2_env TYPE(cp_para_env_type), POINTER :: para_env TYPE(qs_environment_type), POINTER :: qs_env + LOGICAL, INTENT(IN) :: do_kpoints_from_Gamma INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell_3c INTEGER, ALLOCATABLE, DIMENSION(:, :, :), & INTENT(IN) :: cell_to_index_3c @@ -237,9 +242,8 @@ SUBROUTINE compute_mat_P_omega(mat_P_omega, fm_scaled_dm_occ_tau, & fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, & fm_mo_coeff_virt_scaled, mat_dm_occ_global, mat_dm_virt_global, & matrix_s, ispin, & - Eigenval, e_fermi, eps_filter, memory_info, & - unit_nr, & - jquad, do_kpoints_cubic_RPA, qs_env, & + Eigenval, e_fermi, eps_filter, memory_info, unit_nr, & + jquad, do_kpoints_cubic_RPA, do_kpoints_from_Gamma, qs_env, & num_cells_dm, index_to_cell_dm, para_env) ALLOCATE (t_dm_virt(num_cells_dm)) @@ -642,6 +646,7 @@ END SUBROUTINE zero_mat_P_omega !> \param unit_nr ... !> \param jquad ... !> \param do_kpoints_cubic_RPA ... +!> \param do_kpoints_from_Gamma ... !> \param qs_env ... !> \param num_cells_dm ... !> \param index_to_cell_dm ... @@ -651,9 +656,8 @@ SUBROUTINE compute_mat_dm_global(fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, ta fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, & fm_mo_coeff_virt_scaled, mat_dm_occ_global, mat_dm_virt_global, & matrix_s, ispin, & - Eigenval, e_fermi, eps_filter, memory_info, & - unit_nr, & - jquad, do_kpoints_cubic_RPA, qs_env, & + Eigenval, e_fermi, eps_filter, memory_info, unit_nr, & + jquad, do_kpoints_cubic_RPA, do_kpoints_from_Gamma, qs_env, & num_cells_dm, index_to_cell_dm, para_env) TYPE(cp_fm_type), POINTER :: fm_scaled_dm_occ_tau, & @@ -672,7 +676,8 @@ SUBROUTINE compute_mat_dm_global(fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, ta REAL(KIND=dp), INTENT(IN) :: e_fermi, eps_filter LOGICAL, INTENT(IN) :: memory_info INTEGER, INTENT(IN) :: unit_nr, jquad - LOGICAL, INTENT(IN) :: do_kpoints_cubic_RPA + LOGICAL, INTENT(IN) :: do_kpoints_cubic_RPA, & + do_kpoints_from_Gamma TYPE(qs_environment_type), POINTER :: qs_env INTEGER, INTENT(OUT) :: num_cells_dm INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_dm @@ -696,14 +701,25 @@ SUBROUTINE compute_mat_dm_global(fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, ta IF (do_kpoints_cubic_RPA) THEN - CALL compute_transl_dm(mat_dm_occ_global, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, & + CALL compute_transl_dm(mat_dm_occ_global, qs_env, & + ispin, num_integ_points, jquad, e_fermi, tau, & eps_filter, num_cells_dm, index_to_cell_dm, & remove_occ=.FALSE., remove_virt=.TRUE., first_jquad=1) - CALL compute_transl_dm(mat_dm_virt_global, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, & + CALL compute_transl_dm(mat_dm_virt_global, qs_env, & + ispin, num_integ_points, jquad, e_fermi, tau, & eps_filter, num_cells_dm, index_to_cell_dm, & remove_occ=.TRUE., remove_virt=.FALSE., first_jquad=1) + ELSE IF (do_kpoints_from_Gamma) THEN + + CALL compute_mic_dm(mat_dm_occ_global, qs_env, & + ispin, num_integ_points, jquad, e_fermi, tau, & + num_cells_dm, remove_occ=.FALSE., remove_virt=.TRUE., first_jquad=1) + + CALL compute_mic_dm(mat_dm_virt_global, qs_env, & + ispin, num_integ_points, jquad, e_fermi, tau, & + num_cells_dm, remove_occ=.TRUE., remove_virt=.FALSE., first_jquad=1) ELSE num_cells_dm = 1 @@ -824,6 +840,7 @@ SUBROUTINE compute_mat_dm_global(fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, ta END IF ! do kpoints CALL timestop(handle) + END SUBROUTINE compute_mat_dm_global ! ************************************************************************************************** @@ -910,6 +927,7 @@ SUBROUTINE kpoint_density_matrices_rpa(kpoint, tau, e_fermi, remove_occ, remove_ END IF DO i_mo = 1, nmo + IF (ABS(tau*0.5_dp*(eigenvalues(i_mo) - e_fermi)) < stabilize_exp) THEN exp_scaling(i_mo) = EXP(-ABS(tau*(eigenvalues(i_mo) - e_fermi))) ELSE @@ -998,20 +1016,14 @@ SUBROUTINE compute_transl_dm(mat_dm_global, qs_env, ispin, num_integ_points, jqu INTEGER, DIMENSION(3) :: cell_grid_dm TYPE(cell_type), POINTER :: cell TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_global_work, matrix_s_kp - TYPE(dft_control_type), POINTER :: dft_control TYPE(kpoint_type), POINTER :: kpoints TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos - TYPE(neighbor_list_set_p_type), DIMENSION(:), & - POINTER :: sab_nl CALL timeset(routineN, handle) CALL get_qs_env(qs_env, & matrix_s_kp=matrix_s_kp, & - sab_orb=sab_nl, & mos=mos, & - dft_control=dft_control, & - cell=cell, & kpoints=kpoints) nspin = SIZE(mos) @@ -1063,7 +1075,6 @@ SUBROUTINE compute_transl_dm(mat_dm_global, qs_env, ispin, num_integ_points, jqu IF (jquad == first_jquad) THEN NULLIFY (mat_dm_global) -! CALL dbcsr_allocate_matrix_set(mat_dm_global, jquad:num_integ_points, num_cells_dm) ALLOCATE (mat_dm_global(first_jquad:num_integ_points, num_cells_dm)) DO iquad = first_jquad, num_integ_points @@ -1095,6 +1106,227 @@ SUBROUTINE compute_transl_dm(mat_dm_global, qs_env, ispin, num_integ_points, jqu END SUBROUTINE compute_transl_dm +! ************************************************************************************************** +!> \brief ... +!> \param mat_dm_global ... +!> \param qs_env ... +!> \param ispin ... +!> \param num_integ_points ... +!> \param jquad ... +!> \param e_fermi ... +!> \param tau ... +!> \param num_cells_dm ... +!> \param remove_occ ... +!> \param remove_virt ... +!> \param first_jquad ... +! ************************************************************************************************** + SUBROUTINE compute_mic_dm(mat_dm_global, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, & + num_cells_dm, remove_occ, remove_virt, & + first_jquad) + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_global + TYPE(qs_environment_type), POINTER :: qs_env + INTEGER, INTENT(IN) :: ispin, num_integ_points, jquad + REAL(KIND=dp), INTENT(IN) :: e_fermi, tau + INTEGER, INTENT(OUT) :: num_cells_dm + LOGICAL, INTENT(IN) :: remove_occ, remove_virt + INTEGER, INTENT(IN) :: first_jquad + + CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_mic_dm' + + INTEGER :: handle, iquad, jspin, nspin + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_global_work, matrix_s_kp + TYPE(kpoint_type), POINTER :: kpoints_G + TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos + + CALL timeset(routineN, handle) + + NULLIFY (matrix_s_kp, mos) + + CALL get_qs_env(qs_env, & + matrix_s_kp=matrix_s_kp, & + mos=mos) + + kpoints_G => qs_env%mp2_env%ri_rpa_im_time%kpoints_G + + nspin = SIZE(mos) + + num_cells_dm = 1 + + NULLIFY (mat_dm_global_work) + CALL dbcsr_allocate_matrix_set(mat_dm_global_work, nspin, num_cells_dm) + + DO jspin = 1, nspin + + ALLOCATE (mat_dm_global_work(jspin, 1)%matrix) + CALL dbcsr_create(matrix=mat_dm_global_work(jspin, 1)%matrix, & + template=matrix_s_kp(1, 1)%matrix, & + matrix_type=dbcsr_type_no_symmetry) + + CALL dbcsr_reserve_all_blocks(mat_dm_global_work(jspin, 1)%matrix) + + CALL dbcsr_set(mat_dm_global_work(ispin, 1)%matrix, 0.0_dp) + + END DO + + ! density matrices in k-space weighted with EXP(-|e_i-e_F|*t) for occupied orbitals + CALL kpoint_density_matrices_rpa(kpoints_G, tau, e_fermi, & + remove_occ=remove_occ, remove_virt=remove_virt) + + ! density matrices in real space, the cell vectors T for transforming are taken from kpoints%index_to_cell + ! (custom made for RPA) and not from sab_nl (which is symmetric and from SCF) + CALL density_matrix_from_kp_to_mic(kpoints_G, mat_dm_global_work, qs_env) + + ! normally, jquad = 1 to allocate the matrix set, but GW jquad = 0 is the exchange self-energy + IF (jquad == first_jquad) THEN + + NULLIFY (mat_dm_global) + ALLOCATE (mat_dm_global(first_jquad:num_integ_points, num_cells_dm)) + + DO iquad = first_jquad, num_integ_points + NULLIFY (mat_dm_global(iquad, 1)%matrix) + ALLOCATE (mat_dm_global(iquad, 1)%matrix) + CALL dbcsr_create(matrix=mat_dm_global(iquad, 1)%matrix, & + template=matrix_s_kp(1, 1)%matrix, & + matrix_type=dbcsr_type_no_symmetry) + + END DO + + END IF + + CALL dbcsr_copy(mat_dm_global(jquad, 1)%matrix, & + mat_dm_global_work(ispin, 1)%matrix) + + CALL dbcsr_deallocate_matrix_set(mat_dm_global_work) + + CALL timestop(handle) + + END SUBROUTINE compute_mic_dm + + ! ************************************************************************************************** +!> \brief ... +!> \param kpoints_G ... +!> \param mat_dm_global_work ... +!> \param qs_env ... +! ************************************************************************************************** + SUBROUTINE density_matrix_from_kp_to_mic(kpoints_G, mat_dm_global_work, qs_env) + + TYPE(kpoint_type), POINTER :: kpoints_G + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_global_work + TYPE(qs_environment_type), POINTER :: qs_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_from_kp_to_mic' + + INTEGER :: handle, iatom, iatom_old, ik, irow, & + ispin, jatom, jatom_old, jcol, nao, & + ncol_local, nkp, nrow_local, nspin, & + num_cells + INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_from_ao_index + INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices + INTEGER, DIMENSION(:, :), POINTER :: index_to_cell + REAL(KIND=dp) :: contribution, shortest_lattice_vector, & + weight_im, weight_re + REAL(KIND=dp), DIMENSION(2) :: abs_cutoffs_chi_W + REAL(KIND=dp), DIMENSION(3, 3) :: hmat + REAL(KIND=dp), DIMENSION(:), POINTER :: wkp + REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp + TYPE(cell_type), POINTER :: cell + TYPE(cp_fm_type), POINTER :: cpmat, fm_mat_work, rpmat + TYPE(kpoint_env_type), POINTER :: kp + TYPE(mo_set_type), POINTER :: mo_set + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + + CALL timeset(routineN, handle) + + NULLIFY (xkp, wkp) + + NULLIFY (fm_mat_work) + CALL cp_fm_create(fm_mat_work, kpoints_G%kp_env(1)%kpoint_env%wmat(1, 1)%matrix%matrix_struct) + CALL cp_fm_set_all(fm_mat_work, 0.0_dp) + + CALL get_kpoint_info(kpoints_G, nkp=nkp, xkp=xkp, wkp=wkp) + index_to_cell => kpoints_G%index_to_cell + num_cells = SIZE(index_to_cell, 2) + + abs_cutoffs_chi_W = qs_env%mp2_env%ri_rpa_im_time%abs_cutoffs_chi_W + + nspin = SIZE(mat_dm_global_work, 1) + + mo_set => kpoints_G%kp_env(1)%kpoint_env%mos(1, 1)%mo_set + CALL get_mo_set(mo_set, nao=nao) + + ALLOCATE (atom_from_ao_index(nao)) + + CALL get_atom_index_from_basis_function_index(qs_env, atom_from_ao_index, nao, "ORB") + + CALL cp_fm_get_info(matrix=kpoints_G%kp_env(1)%kpoint_env%wmat(1, 1)%matrix, & + nrow_local=nrow_local, & + ncol_local=ncol_local, & + row_indices=row_indices, & + col_indices=col_indices) + + NULLIFY (cell, particle_set) + CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set) + CALL get_cell(cell=cell, h=hmat) + + CALL compute_shortest_lattice_vector(shortest_lattice_vector, hmat, index_to_cell, num_cells) + abs_cutoffs_chi_W(1:2) = qs_env%mp2_env%ri_rpa_im_time%rel_cutoffs_chi_W(1:2)* & + shortest_lattice_vector + + qs_env%mp2_env%ri_rpa_im_time%abs_cutoffs_chi_W(1:2) = abs_cutoffs_chi_W(1:2) + + iatom_old = 0 + jatom_old = 0 + + DO ispin = 1, nspin + + CALL dbcsr_set(mat_dm_global_work(ispin, 1)%matrix, 0.0_dp) + + DO ik = 1, nkp + + kp => kpoints_G%kp_env(ik)%kpoint_env + rpmat => kp%wmat(1, ispin)%matrix + cpmat => kp%wmat(2, ispin)%matrix + + DO irow = 1, nrow_local + DO jcol = 1, ncol_local + + iatom = atom_from_ao_index(row_indices(irow)) + jatom = atom_from_ao_index(col_indices(jcol)) + + IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old) THEN + + CALL compute_weight_re_im(weight_re, weight_im, & + num_cells, iatom, jatom, xkp(1:3, ik), wkp(ik), & + cell, index_to_cell, hmat, particle_set, abs_cutoffs_chi_W) + + iatom_old = iatom + jatom_old = jatom + + END IF + + ! minus sign because of i^2 = -1 + contribution = weight_re*rpmat%local_data(irow, jcol) - & + weight_im*cpmat%local_data(irow, jcol) + + fm_mat_work%local_data(irow, jcol) = fm_mat_work%local_data(irow, jcol) + contribution + + END DO + END DO + + END DO ! ik + + CALL copy_fm_to_dbcsr(fm_mat_work, mat_dm_global_work(ispin, 1)%matrix, keep_sparsity=.FALSE.) + CALL cp_fm_set_all(fm_mat_work, 0.0_dp) + + END DO + + CALL cp_fm_release(fm_mat_work) + DEALLOCATE (atom_from_ao_index) + + CALL timestop(handle) + + END SUBROUTINE density_matrix_from_kp_to_mic + ! ************************************************************************************************** !> \brief ... !> \param kpoints ... @@ -1111,7 +1343,6 @@ SUBROUTINE density_matrix_from_kp_to_transl(kpoints, mat_dm_global_work, index_t INTEGER :: handle, icell, ik, ispin, nkp, nspin, & xcell, ycell, zcell - INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index REAL(KIND=dp) :: arg, coskl, sinkl REAL(KIND=dp), DIMENSION(:), POINTER :: wkp REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp @@ -1121,7 +1352,7 @@ SUBROUTINE density_matrix_from_kp_to_transl(kpoints, mat_dm_global_work, index_t CALL timeset(routineN, handle) - NULLIFY (cell_to_index, xkp, wkp) + NULLIFY (xkp, wkp) NULLIFY (mat_work_re) CALL dbcsr_init_p(mat_work_re) @@ -1135,8 +1366,7 @@ SUBROUTINE density_matrix_from_kp_to_transl(kpoints, mat_dm_global_work, index_t template=mat_dm_global_work(1, 1)%matrix, & matrix_type=dbcsr_type_no_symmetry) - CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, wkp=wkp, & - cell_to_index=cell_to_index) + CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, wkp=wkp) nspin = SIZE(mat_dm_global_work, 1) diff --git a/src/rpa_kpoints.F b/src/rpa_kpoints.F index 309c97b148..abcb997d50 100644 --- a/src/rpa_kpoints.F +++ b/src/rpa_kpoints.F @@ -15,13 +15,17 @@ MODULE rpa_kpoints get_cell,& pbc USE cp_cfm_basic_linalg, ONLY: cp_cfm_cholesky_decompose,& + cp_cfm_column_scale,& cp_cfm_gemm,& - cp_cfm_scale_and_add_fm + cp_cfm_scale_and_add_fm,& + cp_cfm_transpose + USE cp_cfm_diag, ONLY: cp_cfm_heevd 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_type USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm USE cp_fm_types, ONLY: cp_fm_copy_general,& @@ -31,6 +35,7 @@ MODULE rpa_kpoints cp_fm_set_all,& cp_fm_type USE cp_para_types, ONLY: cp_para_env_type + USE cp_units, ONLY: cp_unit_from_cp2k USE dbcsr_api, ONLY: & dbcsr_copy, dbcsr_create, dbcsr_filter, dbcsr_get_block_p, dbcsr_iterator_blocks_left, & dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, & @@ -40,10 +45,12 @@ MODULE rpa_kpoints USE kpoint_types, ONLY: get_kpoint_info,& kpoint_type USE mathconstants, ONLY: gaussi,& + pi,& twopi,& z_one,& z_zero - USE message_passing, ONLY: mp_sum + USE message_passing, ONLY: mp_sum,& + mp_sync USE particle_types, ONLY: particle_type USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type @@ -58,7 +65,7 @@ MODULE rpa_kpoints CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_kpoints' - PUBLIC :: RPA_postprocessing_kp + PUBLIC :: RPA_postprocessing_kp, cp_cfm_robust_cholesky CONTAINS @@ -77,6 +84,7 @@ MODULE rpa_kpoints !> \param wj ... !> \param weights_cos_tf_w_to_t ... !> \param wkp_W ... +!> \param wkp_V ... !> \param do_gw_im_time ... !> \param do_ri_Sigma_x ... !> \param do_kpoints_from_Gamma ... @@ -88,6 +96,7 @@ MODULE rpa_kpoints !> \param mat_P_omega_kp ... !> \param qs_env ... !> \param eps_filter_im_time ... +!> \param unit_nr ... !> \param kpoints ... !> \param fm_mat_L ... !> \param fm_mat_W ... @@ -95,9 +104,10 @@ MODULE rpa_kpoints !> \param mat_SinvVSinv ... ! ************************************************************************************************** SUBROUTINE RPA_postprocessing_kp(dimen_RI, num_integ_points, jquad, nkp, count_ev_sc_GW, para_env, para_env_RPA, & - Erpa, tau_tj, tj, wj, weights_cos_tf_w_to_t, wkp_W, do_gw_im_time, do_ri_Sigma_x, & - do_kpoints_from_Gamma, do_kpoints_cubic_RPA, cfm_mat_W_kp_tau, cfm_mat_Q, ikp_local, & - mat_P_omega, mat_P_omega_kp, qs_env, eps_filter_im_time, kpoints, fm_mat_L, fm_mat_W, & + Erpa, tau_tj, tj, wj, weights_cos_tf_w_to_t, wkp_W, wkp_V, do_gw_im_time, & + do_ri_Sigma_x, do_kpoints_from_Gamma, do_kpoints_cubic_RPA, & + cfm_mat_W_kp_tau, cfm_mat_Q, ikp_local, mat_P_omega, mat_P_omega_kp, & + qs_env, eps_filter_im_time, unit_nr, kpoints, fm_mat_L, fm_mat_W, & fm_mat_RI_global_work, mat_SinvVSinv) INTEGER, INTENT(IN) :: dimen_RI, num_integ_points, jquad, nkp, & @@ -109,7 +119,7 @@ SUBROUTINE RPA_postprocessing_kp(dimen_RI, num_integ_points, jquad, nkp, count_e REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: tj, wj REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & INTENT(IN) :: weights_cos_tf_w_to_t - REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wkp_W + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wkp_W, wkp_V LOGICAL, INTENT(IN) :: do_gw_im_time, do_ri_Sigma_x, & do_kpoints_from_Gamma, & do_kpoints_cubic_RPA @@ -120,6 +130,7 @@ SUBROUTINE RPA_postprocessing_kp(dimen_RI, num_integ_points, jquad, nkp, count_e TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: mat_P_omega, mat_P_omega_kp TYPE(qs_environment_type), POINTER :: qs_env REAL(KIND=dp), INTENT(IN) :: eps_filter_im_time + INTEGER, INTENT(IN) :: unit_nr TYPE(kpoint_type), POINTER :: kpoints TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_mat_L TYPE(cp_fm_p_type), DIMENSION(:), INTENT(IN) :: fm_mat_W @@ -140,7 +151,7 @@ SUBROUTINE RPA_postprocessing_kp(dimen_RI, num_integ_points, jquad, nkp, count_e END IF IF (do_kpoints_from_Gamma) THEN - CALL get_P_cell_T_from_P_gamma(mat_P_omega, qs_env, kpoints, jquad) + CALL get_P_cell_T_from_P_gamma(mat_P_omega, qs_env, kpoints, jquad, unit_nr) END IF CALL transform_P_from_real_space_to_kpoints(mat_P_omega, mat_P_omega_kp, & @@ -148,18 +159,22 @@ SUBROUTINE RPA_postprocessing_kp(dimen_RI, num_integ_points, jquad, nkp, count_e ALLOCATE (trace_Qomega(dimen_RI)) + IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T69,I9,A,I2)') & + 'GW_INFO| Computing chi and W for frequency point: ', jquad, '/', num_integ_points + DO ikp = 1, nkp - ! parallization + ! parallization, we either have all kpoints on all processors or a single kpoint per group IF (ikp_local(ikp) .NE. ikp) CYCLE - ! 1. multiplication Q(iw,k) = K^H(k)P(iw,k)K(k) + ! 1. remove all spurious negative eigenvalues from P(iw,k), multiplication Q(iw,k) = K^H(k)P(iw,k)K(k) CALL compute_Q_kp_RPA(cfm_mat_Q, & mat_P_omega_kp, & fm_mat_L(ikp, 1)%matrix, & fm_mat_L(ikp, 2)%matrix, & fm_mat_RI_global_work, & - dimen_RI, ikp, nkp, ikp_local, para_env) + dimen_RI, ikp, nkp, ikp_local, para_env, & + qs_env%mp2_env%ri_rpa_im_time%make_chi_pos_definite) ! 2. Cholesky decomposition of Id + Q(iw,k) CALL cholesky_decomp_Q(cfm_mat_Q, para_env_RPA, trace_Qomega, dimen_RI) @@ -181,7 +196,7 @@ SUBROUTINE RPA_postprocessing_kp(dimen_RI, num_integ_points, jquad, nkp, count_e fm_mat_L(ikp, 2)%matrix, & dimen_RI, 1, 1, & ikp, tj_dummy, tau_tj_dummy, weights_cos_tf_w_to_t_dummy, & - ikp_local, para_env, kpoints, qs_env, wkp_W, & + ikp_local, para_env, kpoints, qs_env, wkp_V, & mat_SinvVSinv, do_W_and_not_V=.FALSE.) CALL release_dummys(tj_dummy, tau_tj_dummy, weights_cos_tf_w_to_t_dummy) @@ -218,6 +233,63 @@ SUBROUTINE RPA_postprocessing_kp(dimen_RI, num_integ_points, jquad, nkp, count_e END SUBROUTINE RPA_postprocessing_kp +! ************************************************************************************************** +!> \brief Decomposition A = L^H*L of a Hermitian matrix A +!> \param matrix ... +!> \param work ... +!> \param threshold ... +!> \param exponent ... +!> \param min_eigval ... +!> \author Jan Wilhelm, 07.2021 +! ************************************************************************************************** + SUBROUTINE cp_cfm_robust_cholesky(matrix, work, threshold, exponent, min_eigval) + TYPE(cp_cfm_type), POINTER :: matrix, work + REAL(KIND=dp) :: threshold, exponent + REAL(KIND=dp), OPTIONAL :: min_eigval + + CHARACTER(LEN=*), PARAMETER :: routineN = 'cp_cfm_robust_cholesky' + COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), & + czero = CMPLX(0.0_dp, 0.0_dp, KIND=dp) + + COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues_exponent + INTEGER :: handle, i, ncol_global, nrow_global + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues + + CALL timeset(routineN, handle) + + ! Test that matrix is square + CALL cp_cfm_get_info(matrix, nrow_global=nrow_global, ncol_global=ncol_global) + CPASSERT(nrow_global == ncol_global) + ALLOCATE (eigenvalues(nrow_global)) + eigenvalues(:) = 0.0_dp + ALLOCATE (eigenvalues_exponent(nrow_global)) + eigenvalues_exponent(:) = czero + + ! Diagonalize matrix: get eigenvectors and eigenvalues + CALL cp_cfm_heevd(matrix, work, eigenvalues) + + DO i = 1, nrow_global + IF (eigenvalues(i) > threshold) THEN + eigenvalues_exponent(i) = CMPLX((eigenvalues(i))**exponent, threshold, KIND=dp) + ELSE + IF (PRESENT(min_eigval)) THEN + eigenvalues_exponent(i) = CMPLX(min_eigval, 0.0_dp, KIND=dp) + ELSE + eigenvalues_exponent(i) = czero + END IF + END IF + END DO + + CALL cp_cfm_column_scale(work, eigenvalues_exponent) + + CALL cp_cfm_transpose(work, 'C', matrix) + + DEALLOCATE (eigenvalues, eigenvalues_exponent) + + CALL timestop(handle) + + END SUBROUTINE cp_cfm_robust_cholesky + ! ************************************************************************************************** !> \brief ... !> \param cfm_mat_Q ... @@ -230,9 +302,11 @@ END SUBROUTINE RPA_postprocessing_kp !> \param nkp ... !> \param ikp_local ... !> \param para_env ... +!> \param make_chi_pos_definite ... ! ************************************************************************************************** SUBROUTINE compute_Q_kp_RPA(cfm_mat_Q, mat_P_omega_kp, fm_mat_L_re, fm_mat_L_im, & - fm_mat_RI_global_work, dimen_RI, ikp, nkp, ikp_local, para_env) + fm_mat_RI_global_work, dimen_RI, ikp, nkp, ikp_local, para_env, & + make_chi_pos_definite) TYPE(cp_cfm_type), POINTER :: cfm_mat_Q TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: mat_P_omega_kp @@ -241,6 +315,7 @@ SUBROUTINE compute_Q_kp_RPA(cfm_mat_Q, mat_P_omega_kp, fm_mat_L_re, fm_mat_L_im, INTEGER, INTENT(IN) :: dimen_RI, ikp, nkp INTEGER, DIMENSION(:), INTENT(IN) :: ikp_local TYPE(cp_para_env_type), POINTER :: para_env + LOGICAL, INTENT(IN) :: make_chi_pos_definite CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Q_kp_RPA' @@ -262,18 +337,29 @@ SUBROUTINE compute_Q_kp_RPA(cfm_mat_Q, mat_P_omega_kp, fm_mat_L_re, fm_mat_L_im, CALL cp_fm_create(fm_mat_work, fm_mat_L_re%matrix_struct) CALL cp_fm_set_all(fm_mat_work, 0.0_dp) + ! 1. Convert the dbcsr matrix mat_P_omega_kp (that is chi(k,iw)) to a full matrix and + ! distribute it to subgroups CALL mat_P_to_subgroup(mat_P_omega_kp, fm_mat_RI_global_work, & fm_mat_work, cfm_mat_Q, ikp, nkp, ikp_local, para_env) - ! 2. Copy fm_mat_L_re and fm_mat_L_re to cfm_mat_L + ! 2. Remove all negative eigenvalues from chi(k,iw) + IF (make_chi_pos_definite) THEN + CALL cp_cfm_robust_cholesky(cfm_mat_Q, cfm_mat_work, threshold=0.0_dp, exponent=0.5_dp) + CALL cp_cfm_gemm("C", "N", dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_Q, cfm_mat_Q, & + z_zero, cfm_mat_work) + CALL cp_cfm_to_cfm(cfm_mat_work, cfm_mat_Q) + + END IF + + ! 3. Copy fm_mat_L_re and fm_mat_L_re to cfm_mat_L CALL cp_cfm_scale_and_add_fm(z_zero, cfm_mat_L, z_one, fm_mat_L_re) CALL cp_cfm_scale_and_add_fm(z_one, cfm_mat_L, gaussi, fm_mat_L_im) - ! 3. work = P(iw,k)*L(k) + ! 4. work = P(iw,k)*L(k) CALL cp_cfm_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_Q, cfm_mat_L, & z_zero, cfm_mat_work) - ! 4. Q(iw,k) = L^H(k)*work + ! 5. Q(iw,k) = L^H(k)*work CALL cp_cfm_gemm('C', 'N', dimen_RI, dimen_RI, dimen_RI, z_one, cfm_mat_L, cfm_mat_work, & z_zero, cfm_mat_Q) @@ -362,6 +448,8 @@ SUBROUTINE mat_P_to_subgroup(mat_P_omega_kp, fm_mat_RI_global_work, & END IF + CALL mp_sync(para_env%group) + CALL timestop(handle) END SUBROUTINE mat_P_to_subgroup @@ -385,9 +473,18 @@ SUBROUTINE cholesky_decomp_Q(cfm_mat_Q, para_env_RPA, trace_Qomega, dimen_RI) INTEGER :: handle, i_global, iiB, info_chol, & j_global, jjB, ncol_local, nrow_local INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices + TYPE(cp_cfm_type), POINTER :: cfm_mat_Q_tmp, cfm_mat_work CALL timeset(routineN, handle) + NULLIFY (cfm_mat_work) + CALL cp_cfm_create(cfm_mat_work, cfm_mat_Q%matrix_struct) + CALL cp_cfm_set_all(cfm_mat_work, z_zero) + + NULLIFY (cfm_mat_Q_tmp) + CALL cp_cfm_create(cfm_mat_Q_tmp, cfm_mat_Q%matrix_struct) + CALL cp_cfm_set_all(cfm_mat_Q_tmp, z_zero) + ! get info of fm_mat_Q CALL cp_cfm_get_info(matrix=cfm_mat_Q, & nrow_local=nrow_local, & @@ -411,9 +508,15 @@ SUBROUTINE cholesky_decomp_Q(cfm_mat_Q, para_env_RPA, trace_Qomega, dimen_RI) END DO CALL mp_sum(trace_Qomega, para_env_RPA%group) + CALL cp_cfm_to_cfm(cfm_mat_Q, cfm_mat_Q_tmp) + CALL cp_cfm_cholesky_decompose(matrix=cfm_mat_Q, n=dimen_RI, info_out=info_chol) + CPASSERT(info_chol == 0) + CALL cp_cfm_release(cfm_mat_work) + CALL cp_cfm_release(cfm_mat_Q_tmp) + CALL timestop(handle) END SUBROUTINE cholesky_decomp_Q @@ -591,12 +694,13 @@ SUBROUTINE allocate_Wc_kp_tau_GW(cfm_mat_W_kp_tau, cfm_mat_Q, num_integ_points, !> \param qs_env ... !> \param kpoints ... !> \param jquad ... +!> \param unit_nr ... ! ************************************************************************************************** - SUBROUTINE get_P_cell_T_from_P_gamma(mat_P_omega, qs_env, kpoints, jquad) + SUBROUTINE get_P_cell_T_from_P_gamma(mat_P_omega, qs_env, kpoints, jquad, unit_nr) TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN) :: mat_P_omega TYPE(qs_environment_type), POINTER :: qs_env TYPE(kpoint_type), POINTER :: kpoints - INTEGER, INTENT(IN) :: jquad + INTEGER, INTENT(IN) :: jquad, unit_nr CHARACTER(LEN=*), PARAMETER :: routineN = 'get_P_cell_T_from_P_gamma' @@ -604,12 +708,11 @@ SUBROUTINE get_P_cell_T_from_P_gamma(mat_P_omega, qs_env, kpoints, jquad) num_integ_points, row INTEGER, DIMENSION(3) :: cell_grid_P, periodic INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_P - LOGICAL :: found - REAL(KIND=dp) :: cutoff_exp, d_0, sum_exp, weight - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: abs_rab_cell + REAL(KIND=dp) :: abs_rab_cell_i, first_cutoff, & + second_cutoff, weight REAL(KIND=dp), DIMENSION(3) :: cell_vector, rab_cell_i REAL(KIND=dp), DIMENSION(3, 3) :: hmat - REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_to_compute, data_block + REAL(KIND=dp), DIMENSION(:, :), POINTER :: data_block TYPE(cell_type), POINTER :: cell TYPE(dbcsr_iterator_type) :: iter TYPE(particle_type), DIMENSION(:), POINTER :: particle_set @@ -647,59 +750,60 @@ SUBROUTINE get_P_cell_T_from_P_gamma(mat_P_omega, qs_env, kpoints, jquad) mat_P_omega(jquad, 1)%matrix) END DO - ! exponential decay parameter - d_0 = qs_env%mp2_env%ri_rpa_im_time%cutoff - cutoff_exp = 100.0_dp - - ALLOCATE (abs_rab_cell(1:num_cells_P)) + first_cutoff = qs_env%mp2_env%ri_rpa_im_time%abs_cutoffs_chi_W(1) + second_cutoff = qs_env%mp2_env%ri_rpa_im_time%abs_cutoffs_chi_W(2) + + IF (jquad == 1 .AND. unit_nr > 0) THEN + WRITE (unit_nr, '(T3,A,T72,F9.2)') 'GW_INFO| Relative cutoff 1: ', & + qs_env%mp2_env%ri_rpa_im_time%rel_cutoffs_chi_W(1) + 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") + WRITE (unit_nr, '(T3,A,T72,F9.2)') 'GW_INFO| Cutoff radius 2 (Angstroem):', & + 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 + 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' + ELSE + WRITE (unit_nr, '(T3,A,T81)') & + 'GW_INFO| Make chi(iw,k) positive definite? FALSE' + END IF - ! loop over blocks of mat_P_omega(i_cell=1) - CALL dbcsr_iterator_start(iter, mat_P_omega(jquad, 1)%matrix) - DO WHILE (dbcsr_iterator_blocks_left(iter)) - CALL dbcsr_iterator_next_block(iter, row, col, data_block) + END IF - sum_exp = 0.0_dp + DO i_cell = 1, num_cells_P - DO i_cell = 1, num_cells_P + CALL dbcsr_iterator_start(iter, mat_P_omega(jquad, i_cell)%matrix) + DO WHILE (dbcsr_iterator_blocks_left(iter)) + CALL dbcsr_iterator_next_block(iter, row, col, data_block) cell_vector(1:3) = MATMUL(hmat, REAL(index_to_cell_P(1:3, i_cell), dp)) rab_cell_i(1:3) = pbc(particle_set(row)%r(1:3), cell) - & (pbc(particle_set(col)%r(1:3), cell) + cell_vector(1:3)) - abs_rab_cell(i_cell) = SQRT(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2) - IF (abs_rab_cell(i_cell)/d_0 < cutoff_exp) THEN - sum_exp = sum_exp + EXP(-abs_rab_cell(i_cell)/d_0) - END IF - - END DO - - IF (abs_rab_cell(1)/d_0 < cutoff_exp) THEN - weight = EXP(-abs_rab_cell(1)/d_0)/sum_exp - ELSE - weight = 0.0_dp - END IF - data_block(:, :) = data_block(:, :)*weight - - DO i_cell = 2, num_cells_P - - ! THE SYMMETRY EQUIVALENT LATTICE VECTORS ENTER HERE: ONLY LATT. VEC. WITH POS. X-COMP ARE CONSIDERED - IF (abs_rab_cell(i_cell)/d_0 < cutoff_exp .AND. index_to_cell_P(1, i_cell) .LE. 0) THEN - weight = EXP(-abs_rab_cell(i_cell)/d_0)/sum_exp + abs_rab_cell_i = SQRT(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2) + + ! smooth function that is only one if atoms are very close and becomes smoothly zero + ! for larger distances + IF (abs_rab_cell_i < first_cutoff) THEN + weight = 1.0_dp + ELSE IF (abs_rab_cell_i < second_cutoff) THEN + weight = 0.5_dp*(COS(pi*(abs_rab_cell_i - first_cutoff)/ & + (second_cutoff - first_cutoff)) + 1.0_dp) ELSE weight = 0.0_dp END IF - NULLIFY (block_to_compute) - CALL dbcsr_get_block_p(matrix=mat_P_omega(jquad, i_cell)%matrix, & - row=row, col=col, block=block_to_compute, found=found) - CPASSERT(found) - block_to_compute(:, :) = block_to_compute(:, :)*weight + data_block(:, :) = data_block(:, :)*weight END DO + CALL dbcsr_iterator_stop(iter) END DO - CALL dbcsr_iterator_stop(iter) - - DEALLOCATE (abs_rab_cell) CALL timestop(handle) @@ -789,7 +893,6 @@ SUBROUTINE real_space_to_kpoint_transform_rpa(real_mat_kp, imag_mat_kp, mat_real CALL dbcsr_set(real_mat_kp(ik)%matrix, 0.0_dp) CALL dbcsr_set(imag_mat_kp(ik)%matrix, 0.0_dp) - ! JW to check: high memory consumption CALL dbcsr_reserve_all_blocks(real_mat_kp(ik)%matrix) CALL dbcsr_reserve_all_blocks(imag_mat_kp(ik)%matrix) diff --git a/src/rpa_main.F b/src/rpa_main.F index a90d7599ed..01efb470ad 100644 --- a/src/rpa_main.F +++ b/src/rpa_main.F @@ -74,7 +74,11 @@ MODULE rpa_main three_dim_real_array,& two_dim_int_array,& two_dim_real_array - USE qs_environment_types, ONLY: qs_environment_type + USE qs_band_structure, ONLY: calculate_kp_orbitals + USE qs_environment_types, ONLY: get_qs_env,& + qs_env_release,& + qs_environment_type + USE qs_gamma2kp, ONLY: create_kp_from_gamma USE rpa_axk, ONLY: compute_axk_ener USE rpa_gw, ONLY: GW_matrix_operations,& allocate_matrices_gw,& @@ -209,8 +213,8 @@ SUBROUTINE rpa_ri_compute_en(qs_env, Erpa, mp2_env, BIb_C, BIb_C_gw, BIb_C_bse_i INTEGER(KIND=int_8) :: mem INTEGER, ALLOCATABLE, DIMENSION(:) :: dimen_ia, my_ia_end, my_ia_size, & my_ia_start, sub_proc_map, virtual - LOGICAL :: do_minimax_quad, my_open_shell, & - skip_integ_group_opt + LOGICAL :: do_kpoints_from_Gamma, do_minimax_quad, & + my_open_shell, skip_integ_group_opt REAL(KIND=dp) :: allowed_memory, avail_mem, E_Range, Emax, Emin, mem_for_iaK, mem_for_QK, & mem_min, mem_per_group, mem_real, needed_mem REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: BIb_C_2D_bse_ab, BIb_C_2D_bse_ij @@ -288,12 +292,6 @@ SUBROUTINE rpa_ri_compute_en(qs_env, Erpa, mp2_env, BIb_C, BIb_C_gw, BIb_C_bse_i input_integ_group_size = mp2_env%ri_rpa%rpa_integ_group_size IF (my_do_gw .AND. do_minimax_quad) THEN - IF (.NOT. ANY(num_integ_points == (/26, 28, 30, 32, 34/))) THEN - IF (unit_nr > 0) & - CALL cp_warn(__LOCATION__, & - "Please use NUM_INTEG_POINTS 30 for low-scaling GW to "// & - "achieve good accuracy.") - END IF IF (num_integ_points > 34) THEN IF (unit_nr > 0) & CALL cp_warn(__LOCATION__, & @@ -631,6 +629,11 @@ SUBROUTINE rpa_ri_compute_en(qs_env, Erpa, mp2_env, BIb_C, BIb_C_gw, BIb_C_bse_i END IF + do_kpoints_from_Gamma = SUM(mp2_env%ri_rpa_im_time%kp_grid) > 0 + IF (do_kpoints_from_Gamma) THEN + CALL get_bandstruc_and_k_dependent_MOs(qs_env, kpoints, virtual(1)) + END IF + ! Now start the RPA calculation ! fm_mo_coeff_occ, fm_mo_coeff_virt will be deallocated here CALL rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_sub, unit_nr, & @@ -646,7 +649,7 @@ SUBROUTINE rpa_ri_compute_en(qs_env, Erpa, mp2_env, BIb_C, BIb_C_gw, BIb_C_bse_i t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, & starts_array_mc, ends_array_mc, & starts_array_mc_block, ends_array_mc_block, & - matrix_s, kpoints, gd_array, color_sub, & + matrix_s, do_kpoints_from_Gamma, kpoints, gd_array, color_sub, & do_ri_sos_laplace_mp2=do_ri_sos_laplace_mp2) DEALLOCATE (sub_proc_map) @@ -691,6 +694,141 @@ SUBROUTINE rpa_ri_compute_en(qs_env, Erpa, mp2_env, BIb_C, BIb_C_gw, BIb_C_bse_i END SUBROUTINE rpa_ri_compute_en +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param kpoints ... +!> \param virtual ... +! ************************************************************************************************** + SUBROUTINE get_bandstruc_and_k_dependent_MOs(qs_env, kpoints, virtual) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(kpoint_type), POINTER :: kpoints + INTEGER :: virtual + + CHARACTER(LEN=*), PARAMETER :: routineN = 'get_bandstruc_and_k_dependent_MOs' + + INTEGER :: handle, i_dim + INTEGER, DIMENSION(3) :: nkp_grid_G, nkp_grid_W + TYPE(cp_para_env_type), POINTER :: para_env + + CALL timeset(routineN, handle) + + NULLIFY (qs_env%mp2_env%ri_rpa_im_time%kpoints_G, & + qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma, & + qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma_no_xc, & + para_env) + + nkp_grid_W(1:3) = kpoints%nkp_grid(1:3) + DO i_dim = 1, 3 + IF (nkp_grid_W(i_dim) == 1) THEN + nkp_grid_G(i_dim) = 1 + ELSE + nkp_grid_G(i_dim) = nkp_grid_W(i_dim)*qs_env%mp2_env%ri_rpa_im_time%k_mesh_g_factor + END IF + END DO + + 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, & + mp_grid=nkp_grid_G(1:3)) + + CALL timestop(handle) + + END SUBROUTINE get_bandstruc_and_k_dependent_MOs + +! ************************************************************************************************** +!> \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, & + 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, 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' + + INTEGER :: handle + LOGICAL :: my_with_xc_terms + TYPE(qs_environment_type), POINTER :: qs_env_kp_from_Gamma + + 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) + + ! 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) + + CALL qs_env_release(qs_env_kp_from_Gamma) + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param kpgeneral ... +! ************************************************************************************************** + SUBROUTINE get_kpgeneral_for_Sigma_kpoints(qs_env, kpgeneral) + TYPE(qs_environment_type), POINTER :: qs_env + REAL(kind=dp), DIMENSION(:, :), POINTER :: kpgeneral + + CHARACTER(LEN=*), PARAMETER :: routineN = 'get_kpgeneral_for_Sigma_kpoints' + + INTEGER :: handle, i_kp_in_kp_line, i_special_kp, & + ikk, n_kp_in_kp_line, n_special_kp, & + nkp_bandstruc + + CALL timeset(routineN, handle) + + n_special_kp = qs_env%mp2_env%ri_g0w0%n_special_kp + n_kp_in_kp_line = qs_env%mp2_env%ri_g0w0%n_kp_in_kp_line + + nkp_bandstruc = n_kp_in_kp_line*(n_special_kp - 1) + 1 + + ALLOCATE (kpgeneral(3, nkp_bandstruc)) + + kpgeneral(1:3, 1) = qs_env%mp2_env%ri_g0w0%xkp_special_kp(1:3, 1) + + ikk = 1 + + DO i_special_kp = 2, n_special_kp + DO i_kp_in_kp_line = 1, n_kp_in_kp_line + + ikk = ikk + 1 + kpgeneral(1:3, ikk) = qs_env%mp2_env%ri_g0w0%xkp_special_kp(1:3, i_special_kp - 1) + & + REAL(i_kp_in_kp_line, KIND=dp)/REAL(n_kp_in_kp_line, KIND=dp)* & + (qs_env%mp2_env%ri_g0w0%xkp_special_kp(1:3, i_special_kp) - & + qs_env%mp2_env%ri_g0w0%xkp_special_kp(1:3, i_special_kp - 1)) + + END DO + END DO + + CALL timestop(handle) + + END SUBROUTINE get_kpgeneral_for_Sigma_kpoints + ! ************************************************************************************************** !> \brief reorder the local data in such a way to help the next stage of matrix creation; !> now the data inside the group are divided into a ia x K matrix (BIb_C_2D); @@ -1158,6 +1296,7 @@ END SUBROUTINE create_occ_virt_mo_coeffs !> \param starts_array_mc_block ... !> \param ends_array_mc_block ... !> \param matrix_s ... +!> \param do_kpoints_from_Gamma ... !> \param kpoints ... !> \param gd_array ... !> \param color_sub ... @@ -1176,7 +1315,7 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, & starts_array_mc, ends_array_mc, & starts_array_mc_block, ends_array_mc_block, & - matrix_s, kpoints, gd_array, color_sub, & + matrix_s, do_kpoints_from_Gamma, kpoints, gd_array, color_sub, & do_ri_sos_laplace_mp2) TYPE(qs_environment_type), POINTER :: qs_env @@ -1216,6 +1355,7 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s starts_array_mc_block, & ends_array_mc_block TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s + LOGICAL :: do_kpoints_from_Gamma TYPE(kpoint_type), POINTER :: kpoints TYPE(group_dist_d1_type), INTENT(IN) :: gd_array INTEGER, INTENT(IN) :: color_sub @@ -1236,15 +1376,15 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s INTEGER, DIMENSION(:), POINTER :: col_blk_size, prim_blk_sizes, & RI_blk_sizes LOGICAL :: do_apply_ic_corr_to_gw, do_gw_im_time, do_ic_model, do_kpoints_cubic_RPA, & - do_kpoints_from_Gamma, do_periodic, do_ri_Sigma_x, exit_ev_gw, first_cycle, & - first_cycle_periodic_correction, my_open_shell, print_ic_values + do_periodic, do_ri_Sigma_x, exit_ev_gw, first_cycle, first_cycle_periodic_correction, & + my_open_shell, print_ic_values LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :) :: has_mat_P_blocks REAL(KIND=dp) :: a_scaling, alpha, dbcsr_time, e_axk, e_axk_corr, eps_filter, & eps_filter_im_time, eps_min_trans, ext_scaling, fermi_level_offset, & fermi_level_offset_input, my_flop_rate, omega, omega_max_fit, omega_old, tau, tau_old REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: delta_corr, e_fermi, tau_tj, tau_wj, tj, & trace_Qomega, vec_omega_fit_gw, wj, & - wkp_W + wkp_V, wkp_W REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Eigenval_kp, Eigenval_last, Eigenval_scf, & Eigenval_scf_kp, vec_W_gw, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, & weights_sin_tf_t_to_w @@ -1288,7 +1428,6 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s print_ic_values = mp2_env%ri_g0w0%print_ic_values do_periodic = mp2_env%ri_g0w0%do_periodic do_kpoints_cubic_RPA = mp2_env%ri_rpa_im_time%do_im_time_kpoints - do_kpoints_from_Gamma = SUM(mp2_env%ri_rpa_im_time%kp_grid) > 0 ! For SOS-MP2 only gemm is implemented mm_style = wfc_mm_style_gemm @@ -1319,16 +1458,6 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s IF (do_im_time) THEN CPASSERT(do_minimax_quad .OR. do_ri_sos_laplace_mp2) - END IF - - IF (do_ic_model) THEN - ! image charge model only implemented for cubic scaling GW - CPASSERT(do_gw_im_time) - CPASSERT(.NOT. do_periodic) - END IF - - ! set up the least-square time grid and other matrices specifically for imag time - IF (do_im_time) THEN group_size_P = mp2_env%ri_rpa_im_time%group_size_P cut_memory = mp2_env%ri_rpa_im_time%cut_memory @@ -1349,7 +1478,7 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s col_blk_size, & do_ic_model, do_kpoints_cubic_RPA, & do_kpoints_from_Gamma, do_ri_Sigma_x, my_open_shell, & - has_mat_P_blocks, wkp_W, & + has_mat_P_blocks, wkp_W, wkp_V, & cfm_mat_Q, fm_mat_L, fm_mat_RI_global_work, fm_mat_work, fm_mo_coeff_occ_scaled, & fm_mo_coeff_virt_scaled, mat_dm, mat_L, mat_M_P_munu_occ, mat_M_P_munu_virt, mat_P_global_copy, & mat_SinvVSinv, mat_P_omega, mat_P_omega_kp, mat_work) @@ -1383,6 +1512,13 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s END IF + IF (do_ic_model) THEN + ! image charge model only implemented for cubic scaling GW + CPASSERT(do_gw_im_time) + CPASSERT(.NOT. do_periodic) + IF (cut_memory .NE. 1) CPABORT("For IC, use MEMORY_CUT 1 in the LOW_SCALING section.") + END IF + ALLOCATE (e_fermi(nspins)) IF (do_minimax_quad .OR. do_ri_sos_laplace_mp2) THEN CALL get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_integ_points, do_im_time, & @@ -1499,7 +1635,8 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s eps_filter_im_time, Eigenval(:, 1), nmo, & num_integ_points, cut_memory, & unit_nr, mp2_env, para_env, & - qs_env, index_to_cell_3c, cell_to_index_3c, & + qs_env, do_kpoints_from_Gamma, & + index_to_cell_3c, cell_to_index_3c, & has_mat_P_blocks, do_ri_sos_laplace_mp2, & dbcsr_time, dbcsr_nflop) @@ -1522,7 +1659,8 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s eps_filter_im_time, Eigenval(:, 2), nmo, & num_integ_points, cut_memory, & unit_nr, mp2_env, para_env, & - qs_env, index_to_cell_3c, cell_to_index_3c, & + qs_env, do_kpoints_from_Gamma, & + index_to_cell_3c, cell_to_index_3c, & has_mat_P_blocks, do_ri_sos_laplace_mp2, & dbcsr_time, dbcsr_nflop) ELSE @@ -1539,7 +1677,8 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s eps_filter_im_time, Eigenval(:, 2), nmo, & num_integ_points, cut_memory, & unit_nr, mp2_env, para_env, & - qs_env, index_to_cell_3c, cell_to_index_3c, & + qs_env, do_kpoints_from_Gamma, & + index_to_cell_3c, cell_to_index_3c, & has_mat_P_blocks, do_ri_sos_laplace_mp2, & dbcsr_time, dbcsr_nflop) END IF ! do_ri_sos_laplace_mp2 @@ -1619,10 +1758,10 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN CALL RPA_postprocessing_kp(dimen_RI, num_integ_points, jquad, nkp, count_ev_sc_GW, para_env, & - para_env_RPA, Erpa, tau_tj, tj, wj, weights_cos_tf_w_to_t, wkp_W, do_gw_im_time, & - do_ri_Sigma_x, do_kpoints_from_Gamma, do_kpoints_cubic_RPA, cfm_mat_W_kp_tau, & - cfm_mat_Q, ikp_local, mat_P_omega(:, :, 1), & - mat_P_omega_kp, qs_env, eps_filter_im_time, & + para_env_RPA, Erpa, tau_tj, tj, wj, weights_cos_tf_w_to_t, & + wkp_W, wkp_V, do_gw_im_time, do_ri_Sigma_x, do_kpoints_from_Gamma, & + do_kpoints_cubic_RPA, cfm_mat_W_kp_tau, cfm_mat_Q, ikp_local, & + mat_P_omega(:, :, 1), mat_P_omega_kp, qs_env, eps_filter_im_time, unit_nr, & kpoints, fm_mat_L, fm_mat_W, fm_mat_RI_global_work, mat_SinvVSinv) ELSE CALL RPA_postprocessing_nokp(dimen_RI_red, trace_Qomega, fm_mat_Q(1)%matrix, para_env_RPA, Erpa, wj(jquad)) @@ -1791,11 +1930,13 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s cell_to_index_3c, do_ic_model, & do_kpoints_cubic_RPA, do_kpoints_from_Gamma, do_ri_Sigma_x, & has_mat_P_blocks, & - wkp_W, cfm_mat_Q, fm_mat_L, fm_mat_RI_global_work, fm_mat_work, & + wkp_W, wkp_V, cfm_mat_Q, fm_mat_L, fm_mat_RI_global_work, fm_mat_work, & fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm, mat_L, & mat_SinvVSinv, mat_P_omega, mat_P_omega_kp, & - t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, & - mat_work) + t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, mat_work, & + qs_env%mp2_env%ri_rpa_im_time%kpoints_G, & + qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma, & + qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma_no_xc) IF (my_do_gw) THEN CALL deallocate_matrices_gw_im_time(weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, do_ic_model, & diff --git a/src/rpa_util.F b/src/rpa_util.F index 0acec11321..c5cd82bf83 100644 --- a/src/rpa_util.F +++ b/src/rpa_util.F @@ -49,6 +49,7 @@ MODULE rpa_util wfc_mm_style_syrk USE kinds, ONLY: dp USE kpoint_types, ONLY: get_kpoint_info,& + kpoint_release,& kpoint_type USE mathconstants, ONLY: z_zero USE message_passing, ONLY: mp_sum @@ -105,6 +106,7 @@ MODULE rpa_util !> \param my_open_shell ... !> \param has_mat_P_blocks ... !> \param wkp_W ... +!> \param wkp_V ... !> \param cfm_mat_Q ... !> \param fm_mat_L ... !> \param fm_mat_RI_global_work ... @@ -132,7 +134,7 @@ SUBROUTINE alloc_im_time(qs_env, para_env, dimen_RI, dimen_RI_red, num_integ_poi col_blk_size, & do_ic_model, do_kpoints_cubic_RPA, & do_kpoints_from_Gamma, do_ri_Sigma_x, my_open_shell, & - has_mat_P_blocks, wkp_W, & + has_mat_P_blocks, wkp_W, wkp_V, & cfm_mat_Q, fm_mat_L, fm_mat_RI_global_work, fm_mat_work, fm_mo_coeff_occ_scaled, & fm_mo_coeff_virt_scaled, mat_dm, mat_L, mat_M_P_munu_occ, mat_M_P_munu_virt, mat_P_global_copy, & mat_SinvVSinv, mat_P_omega, mat_P_omega_kp, & @@ -164,7 +166,7 @@ SUBROUTINE alloc_im_time(qs_env, para_env, dimen_RI, dimen_RI_red, num_integ_poi LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :), & INTENT(OUT) :: has_mat_P_blocks REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & - INTENT(OUT) :: wkp_W + INTENT(OUT) :: wkp_W, wkp_V TYPE(cp_cfm_type), POINTER :: cfm_mat_Q TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_mat_L TYPE(cp_fm_type), POINTER :: fm_mat_RI_global_work, fm_mat_work, & @@ -217,7 +219,7 @@ SUBROUTINE alloc_im_time(qs_env, para_env, dimen_RI, dimen_RI_red, num_integ_poi IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN CALL get_sub_para_kp(fm_struct_sub_kp, para_env, kpoints%nkp, & - dimen_RI, ikp_local, first_ikp_local, do_kpoints_cubic_RPA) + dimen_RI, ikp_local, first_ikp_local) NULLIFY (cfm_mat_Q) CALL cp_cfm_create(cfm_mat_Q, fm_struct_sub_kp) @@ -237,8 +239,8 @@ SUBROUTINE alloc_im_time(qs_env, para_env, dimen_RI, dimen_RI_red, num_integ_poi CALL get_kpoint_info(kpoints, nkp=nkp) ! compute k-point weights such that functions 1/k^2, 1/k and const function are ! integrated correctly - CALL compute_wkp_W(wkp_W, kpoints, cell%hmat, cell%h_inv, & - qs_env%mp2_env%ri_rpa_im_time%exp_kpoints, periodic) + CALL compute_wkp_W(qs_env, wkp_W, wkp_V, kpoints, cell%h_inv, periodic) + ELSE nkp = 1 END IF @@ -914,16 +916,14 @@ END SUBROUTINE RPA_postprocessing_nokp !> \param dimen_RI ... !> \param ikp_local ... !> \param first_ikp_local ... -!> \param do_kpoints_cubic_RPA ... ! ************************************************************************************************** SUBROUTINE get_sub_para_kp(fm_struct_sub_kp, para_env, nkp, dimen_RI, & - ikp_local, first_ikp_local, do_kpoints_cubic_RPA) + ikp_local, first_ikp_local) TYPE(cp_fm_struct_type), POINTER :: fm_struct_sub_kp TYPE(cp_para_env_type), POINTER :: para_env INTEGER, INTENT(IN) :: nkp, dimen_RI INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: ikp_local INTEGER, INTENT(OUT) :: first_ikp_local - LOGICAL, INTENT(IN) :: do_kpoints_cubic_RPA CHARACTER(len=*), PARAMETER :: routineN = 'get_sub_para_kp' @@ -934,13 +934,9 @@ SUBROUTINE get_sub_para_kp(fm_struct_sub_kp, para_env, nkp, dimen_RI, & CALL timeset(routineN, handle) - IF (nkp > para_env%num_pe .OR. do_kpoints_cubic_RPA) THEN - ! we have all kpoints on every processpr - num_proc_per_kp = para_env%num_pe - ELSE - ! we have only one kpoint per group - num_proc_per_kp = para_env%num_pe/nkp - END IF + ! JW we use all processors for every k-point, a subgroup does not improve the + ! performance substantially + num_proc_per_kp = para_env%num_pe color_sub_kp = MOD(para_env%mepos/num_proc_per_kp, nkp) CALL cp_para_env_split(para_env_sub_kp, para_env, color_sub_kp) @@ -960,10 +956,7 @@ SUBROUTINE get_sub_para_kp(fm_struct_sub_kp, para_env, nkp, dimen_RI, & ikp_local = 0 first_ikp_local = 1 DO ikp = 1, nkp - IF (nkp > para_env%num_pe .OR. do_kpoints_cubic_RPA .OR. ikp == color_sub_kp + 1) THEN - ikp_local(ikp) = ikp - first_ikp_local = ikp - END IF + ikp_local(ikp) = ikp END DO CALL timestop(handle) @@ -985,6 +978,7 @@ END SUBROUTINE get_sub_para_kp !> \param do_ri_Sigma_x ... !> \param has_mat_P_blocks ... !> \param wkp_W ... +!> \param wkp_V ... !> \param cfm_mat_Q ... !> \param fm_mat_L ... !> \param fm_mat_RI_global_work ... @@ -1001,16 +995,20 @@ END SUBROUTINE get_sub_para_kp !> \param t_3c_O_compressed ... !> \param t_3c_O_ind ... !> \param mat_work ... +!> \param kpoints_G ... +!> \param kpoints_Sigma ... +!> \param kpoints_Sigma_no_xc ... ! ************************************************************************************************** SUBROUTINE dealloc_im_time(fm_mo_coeff_occ, fm_mo_coeff_virt, & fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, ikp_local, index_to_cell_3c, & cell_to_index_3c, do_ic_model, & do_kpoints_cubic_RPA, do_kpoints_from_Gamma, do_ri_Sigma_x, & has_mat_P_blocks, & - wkp_W, cfm_mat_Q, fm_mat_L, fm_mat_RI_global_work, fm_mat_work, & + wkp_W, wkp_V, cfm_mat_Q, fm_mat_L, fm_mat_RI_global_work, fm_mat_work, & fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm, mat_L, & mat_SinvVSinv, mat_P_omega, mat_P_omega_kp, & - t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, mat_work) + t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, & + mat_work, kpoints_G, kpoints_Sigma, kpoints_Sigma_no_xc) TYPE(cp_fm_p_type), DIMENSION(:) :: fm_mo_coeff_occ, fm_mo_coeff_virt TYPE(cp_fm_type), POINTER :: fm_scaled_dm_occ_tau, & @@ -1025,7 +1023,7 @@ SUBROUTINE dealloc_im_time(fm_mo_coeff_occ, fm_mo_coeff_virt, & LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :), & INTENT(INOUT) :: has_mat_P_blocks REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & - INTENT(INOUT) :: wkp_W + INTENT(INOUT) :: wkp_W, wkp_V TYPE(cp_cfm_type), POINTER :: cfm_mat_Q TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: fm_mat_L TYPE(cp_fm_type), POINTER :: fm_mat_RI_global_work, fm_mat_work, & @@ -1042,6 +1040,8 @@ SUBROUTINE dealloc_im_time(fm_mo_coeff_occ, fm_mo_coeff_virt, & TYPE(block_ind_type), ALLOCATABLE, & DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_O_ind TYPE(dbcsr_type), POINTER :: mat_work + TYPE(kpoint_type), POINTER :: kpoints_G, kpoints_Sigma, & + kpoints_Sigma_no_xc CHARACTER(LEN=*), PARAMETER :: routineN = 'dealloc_im_time' @@ -1120,6 +1120,7 @@ SUBROUTINE dealloc_im_time(fm_mo_coeff_occ, fm_mo_coeff_virt, & CALL cp_fm_release(fm_mat_RI_global_work) CALL dbcsr_deallocate_matrix_set(mat_P_omega_kp) DEALLOCATE (wkp_W) + DEALLOCATE (wkp_V) END IF cut_memory = SIZE(t_3c_O_compressed, 3) @@ -1134,6 +1135,15 @@ SUBROUTINE dealloc_im_time(fm_mo_coeff_occ, fm_mo_coeff_virt, & END DO DEALLOCATE (t_3c_O_compressed) + IF (do_kpoints_from_Gamma) THEN + IF (ASSOCIATED(kpoints_G)) CALL kpoint_release(kpoints_G) + ! JW later for band structures from GW + IF (.NOT. do_kpoints_from_Gamma) THEN + IF (ASSOCIATED(kpoints_Sigma)) CALL kpoint_release(kpoints_Sigma) + IF (ASSOCIATED(kpoints_Sigma_no_xc)) CALL kpoint_release(kpoints_Sigma_no_xc) + END IF + END IF + CALL timestop(handle) END SUBROUTINE dealloc_im_time diff --git a/tests/QS/regtest-gw-cubic/G0W0_H2O_PBE_periodic.inp b/tests/QS/regtest-gw-cubic/G0W0_H2O_PBE_periodic.inp index 5b120d0087..0f034af7a2 100644 --- a/tests/QS/regtest-gw-cubic/G0W0_H2O_PBE_periodic.inp +++ b/tests/QS/regtest-gw-cubic/G0W0_H2O_PBE_periodic.inp @@ -60,11 +60,11 @@ ANALYTIC_CONTINUATION TWO_POLE CROSSING_SEARCH Z_SHOT RI_SIGMA_X FALSE - PERIODIC - &PERIODIC + PERIODIC_CORRECTION + &PERIODIC_CORRECTION # That value should be chosen much larger ! NUM_OMEGA_POINTS 5 - &END + &END PERIODIC_CORRECTION &END GW &END RI_RPA MEMORY 200. diff --git a/tests/QS/regtest-gw-cubic/G0W0_kpoints_from_Gamma_RI_regularization_1E-3.inp b/tests/QS/regtest-gw-cubic/G0W0_kpoints_from_Gamma_RI_regularization_1E-3.inp new file mode 100644 index 0000000000..8901c9d7b2 --- /dev/null +++ b/tests/QS/regtest-gw-cubic/G0W0_kpoints_from_Gamma_RI_regularization_1E-3.inp @@ -0,0 +1,83 @@ +&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 10000000 + &END SCF + &XC + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &WF_CORRELATION + &INTEGRALS + SIZE_LATTICE_SUM 5 + &END INTEGRALS + &LOW_SCALING + KPOINTS 1 1 4 + REGULARIZATION_RI 1.0E-3 + &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 + &END GW + &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 Z + &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-cubic/TEST_FILES b/tests/QS/regtest-gw-cubic/TEST_FILES index 06e101f668..f6bcd0a92b 100644 --- a/tests/QS/regtest-gw-cubic/TEST_FILES +++ b/tests/QS/regtest-gw-cubic/TEST_FILES @@ -7,6 +7,7 @@ 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.47 G0W0_OH_PBE.inp 79 1e-05 11.80 -G0W0_kpoints_from_Gamma.inp 78 1e-05 15.20 +G0W0_kpoints_from_Gamma.inp 78 1e-05 14.55 +G0W0_kpoints_from_Gamma_RI_regularization_1E-3.inp 78 1e-05 14.71 G0W0_OH_PBE_svd.inp 79 1e-05 11.80 #EOF diff --git a/tests/QS/regtest-gw/G0W0_H2O_PBE_periodic.inp b/tests/QS/regtest-gw/G0W0_H2O_PBE_periodic.inp index 2ddef3e1f8..8b6e7d7e20 100644 --- a/tests/QS/regtest-gw/G0W0_H2O_PBE_periodic.inp +++ b/tests/QS/regtest-gw/G0W0_H2O_PBE_periodic.inp @@ -66,13 +66,14 @@ EV_GW_ITER 3 UPDATE_XC_ENERGY RI_SIGMA_X FALSE - PERIODIC - &PERIODIC + ! the periodic correction is an old method that is not recommended to use + PERIODIC_CORRECTION + &PERIODIC_CORRECTION KPOINTS 12 12 12 NUM_KP_GRIDS 1 EXTRAPOLATE DO_AUX_BAS_GW - &END PERIODIC + &END PERIODIC_CORRECTION &END GW &END RI_RPA MEMORY 200.