diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 2d44f742732..9dbc1e17520 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -889,6 +889,7 @@ OBJS_LCAO=evolve_elec.o\ center2orb_orb22.o\ OBJS_MODULE_RI=conv_coulomb_pot_k.o\ + smooth_cut_coulomb.o\ invmat_ri.o\ lri_cv.o\ lri_cv_tools.o\ diff --git a/source/source_hamilt/module_xc/exx_info.cpp b/source/source_hamilt/module_xc/exx_info.cpp index 98896c44900..56da05e41be 100644 --- a/source/source_hamilt/module_xc/exx_info.cpp +++ b/source/source_hamilt/module_xc/exx_info.cpp @@ -8,6 +8,8 @@ #include #include #include +#include +#include #include //---------------------------------------------------------- @@ -244,6 +246,32 @@ bool init_exx_info(Exx_Info& exx_info, const Input_para& inp) // nspin=4/SOC via restore_dm + restore_HR_nspin4), so symmetry=1 is honored here. exx_info.sync_from_global(); + // Explicit INPUT -> RI configuration; only init_spencer consumes these + // tags, so the full-Coulomb/ground-state exchange kernel is unchanged. + if (inp.rpa) + { + const bool cut_nondefault = inp.rpa_cut_coulomb != "hard" || inp.rpa_cut_radius != "legacy"; + if (cut_nondefault + && exx_info.info_ri.coulomb_param.count(Conv_Coulomb_Pot_K::Coulomb_Type::Fock) == 0) + { + throw std::invalid_argument( + "Non-default rpa_cut_coulomb/rpa_cut_radius requires a Fock Coulomb component; " + "Erfc-only configurations such as default HSE are unsupported."); + } + for (auto& cut_family : exx_info.info_ri.coulomb_param) + { + if (cut_family.first != Conv_Coulomb_Pot_K::Coulomb_Type::Fock) { continue; } + for (auto& cut_parameter : cut_family.second) + { + cut_parameter["rpa_cut_coulomb"] = inp.rpa_cut_coulomb; + cut_parameter["rpa_cut_radius"] = inp.rpa_cut_radius; + std::ostringstream cut_width_text; + cut_width_text.precision(17); + cut_width_text << inp.rpa_cut_width; + cut_parameter["rpa_cut_width"] = cut_width_text.str(); + } + } + } } return generate_opt_orb; diff --git a/source/source_hamilt/module_xc/test/CMakeLists.txt b/source/source_hamilt/module_xc/test/CMakeLists.txt index 82e8188923c..4a51c38880d 100644 --- a/source/source_hamilt/module_xc/test/CMakeLists.txt +++ b/source/source_hamilt/module_xc/test/CMakeLists.txt @@ -1,3 +1,9 @@ +AddTest( + TARGET MODULE_HAMILT_RpaCutoffConfig + LIBS parameter base device + SOURCES test_exx_info.cpp ../exx_info.cpp +) + AddTest( TARGET MODULE_HAMILT_XCTest_PBE LIBS parameter MPI::MPI_CXX Libxc::xc diff --git a/source/source_hamilt/module_xc/test/test_exx_info.cpp b/source/source_hamilt/module_xc/test/test_exx_info.cpp new file mode 100644 index 00000000000..8788d56240d --- /dev/null +++ b/source/source_hamilt/module_xc/test/test_exx_info.cpp @@ -0,0 +1,61 @@ +#include "../exx_info.h" +#include "source_io/module_parameter/input_parameter.h" +#include +#include + +TEST(RpaCutoffConfig, HseDefaultsRemainSupportedButNondefaultCutoffsAreRejected) +{ + Input_para cut_input; + cut_input.basis_type = "lcao"; + cut_input.rpa = true; + cut_input.dft_functional = "hse"; + // Prepared HSE defaults, after INPUT reset hooks have resolved "default". + cut_input.exx_fock_alpha.clear(); + cut_input.exx_fock_lambda.clear(); + cut_input.exx_erfc_alpha = {"0.25"}; + cut_input.exx_erfc_omega = {"0.11"}; + cut_input.exx_singularity_correction = "limits"; + cut_input.exx_real_number = "1"; + cut_input.exx_ccp_rmesh_times = "2"; + Exx_Info cut_default_info; + EXPECT_NO_THROW(init_exx_info(cut_default_info, cut_input)); + EXPECT_EQ(cut_default_info.info_ri.coulomb_param.count(Conv_Coulomb_Pot_K::Coulomb_Type::Fock), 0); + EXPECT_EQ(cut_default_info.info_ri.coulomb_param.count(Conv_Coulomb_Pot_K::Coulomb_Type::Erfc), 1); + + cut_input.rpa_cut_coulomb = "log_erfc"; + Exx_Info cut_smooth_info; + EXPECT_THROW(init_exx_info(cut_smooth_info, cut_input), std::invalid_argument); + cut_input.rpa_cut_coulomb = "hard"; + cut_input.rpa_cut_radius = "inscribed"; + Exx_Info cut_inscribed_info; + EXPECT_THROW(init_exx_info(cut_inscribed_info, cut_input), std::invalid_argument); +} + +TEST(RpaCutoffConfig, PbeRiInitializationMatchesLibriAvailability) +{ + Input_para cut_input; + cut_input.basis_type = "lcao"; + cut_input.rpa = true; + cut_input.dft_functional = "pbe"; + cut_input.exx_fock_lambda.clear(); + cut_input.exx_erfc_omega.clear(); + cut_input.exx_singularity_correction = "limits"; + cut_input.exx_real_number = "1"; + cut_input.exx_ccp_rmesh_times = "2"; + cut_input.rpa_cut_coulomb = "log_erfc"; + cut_input.rpa_cut_radius = "inscribed"; + Exx_Info cut_pbe_info; + ASSERT_NO_THROW(init_exx_info(cut_pbe_info, cut_input)); +#ifdef __EXX + // ENABLE_LIBRI defines __EXX and enables PBE RPA initialization of info_ri. + ASSERT_EQ(cut_pbe_info.info_ri.coulomb_param.count(Conv_Coulomb_Pot_K::Coulomb_Type::Fock), 1); + const auto& cut_fock = cut_pbe_info.info_ri.coulomb_param.at(Conv_Coulomb_Pot_K::Coulomb_Type::Fock); + ASSERT_EQ(cut_fock.size(), 1); + EXPECT_EQ(cut_fock.front().at("rpa_cut_coulomb"), "log_erfc"); + EXPECT_EQ(cut_fock.front().at("rpa_cut_radius"), "inscribed"); + EXPECT_DOUBLE_EQ(std::stod(cut_fock.front().at("rpa_cut_width")), cut_input.rpa_cut_width); +#else + // Keep coverage when LibRI is disabled: PBE must skip RI initialization. + EXPECT_TRUE(cut_pbe_info.info_ri.coulomb_param.empty()); +#endif +} diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 1caf2ec975e..3df93c47515 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -677,6 +677,9 @@ struct Input_para bool exx_symmetry_realspace = true; ///< whether to reduce the real-space sector in when using symmetry=1 in EXX calculation double rpa_ccp_rmesh_times = 10.0; ///< how many times larger the radial mesh required for ///< calculating Columb potential is to that of atomic orbitals + std::string rpa_cut_coulomb = "hard"; ///< post-SCF cut kernel: hard or log_erfc + std::string rpa_cut_radius = "legacy"; ///< legacy radius or BvK inscribed radius + double rpa_cut_width = 1.091868692359; ///< multiplicative log-erfc width double exx_cs_inv_thr = -1; ///< threshold to inverse Vq in abfs for generating Cs double shrink_abfs_pca_thr = -1; ///< threshold to shrink auxiliary basis for GW/RPA double shrink_LU_inv_thr = 1e-6; ///< threshold to get inverse of overlap matrix by LU decomposition in auxiliary basis representation diff --git a/source/source_io/module_parameter/read_inp_exx_dftu.cpp b/source/source_io/module_parameter/read_inp_exx_dftu.cpp index 9428c0e7a40..187d710e3cc 100644 --- a/source/source_io/module_parameter/read_inp_exx_dftu.cpp +++ b/source/source_io/module_parameter/read_inp_exx_dftu.cpp @@ -569,6 +569,42 @@ Use the same scheme for consistent SCF and NSCF bands. A source configuration mi }; this->add_item(item); } + { + Input_Item item("rpa_cut_coulomb"); + item.annotation = "post-SCF truncated Coulomb kernel for RPA/GW export"; + item.category = "Exact Exchange (LCAO)"; + item.type = "String"; + item.description = "Post-SCF cut-Coulomb kernel: hard preserves Spencer truncation; log_erfc uses erfc(log(r/Rc)/log(rpa_cut_width))/2r. Applies to Fock components in compressed and uncompressed cut exports; Erfc components retain their existing kernel. Non-default kernel/radius requests require a Fock component. Full Coulomb is unchanged. Forces and stresses for log_erfc are not validated."; + item.default_value = "hard"; + item.unit = ""; + read_sync_string(input.rpa_cut_coulomb); + item.domain.set_choices({"hard", "log_erfc"}); + this->add_item(item); + } + { + Input_Item item("rpa_cut_radius"); + item.annotation = "radius convention for post-SCF cut Coulomb"; + item.category = "Exact Exchange (LCAO)"; + item.type = "String"; + item.description = "legacy preserves the existing cutoff radius. inscribed uses half the shortest nonzero BvK superlattice vector, including the lattice scale in Bohr. Applies only to Fock components; inscribed is rejected if the RPA Coulomb configuration has no Fock component. Independent of rpa_cut_coulomb; intended for controlled radius-versus-smoothing comparisons."; + item.default_value = "legacy"; + item.unit = ""; + read_sync_string(input.rpa_cut_radius); + item.domain.set_choices({"legacy", "inscribed"}); + this->add_item(item); + } + { + Input_Item item("rpa_cut_width"); + item.annotation = "multiplicative width of the log-erfc cutoff"; + item.category = "Exact Exchange (LCAO)"; + item.type = "Real"; + item.description = "Width factor for log_erfc. Must be finite and greater than 1 and at most 1.25. The default matches the logged Si FHI-aims width; it is not a universal FHI-aims constant. The finite radial integration range must resolve the smooth tail."; + item.default_value = "1.091868692359"; + item.unit = ""; + read_sync_double(input.rpa_cut_width); + item.domain.set_rule("value>1 and value<=1.25"); + this->add_item(item); + } { Input_Item item("exx_symmetry_realspace"); item.annotation = "whether to reduce real-space sector in Hexx calculation"; diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp index b1221ff0b77..eb2f04d5594 100644 --- a/source/source_io/test_serial/read_input_item_test.cpp +++ b/source/source_io/test_serial/read_input_item_test.cpp @@ -2072,6 +2072,30 @@ TEST_F(InputTest, Item_test2) output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output, testing::HasSubstr("NOTICE")); } + { // independent smoothing and radius switches + const auto cut_mode_item = find_label("rpa_cut_coulomb", readinput.input_lists); + const auto cut_radius_item = find_label("rpa_cut_radius", readinput.input_lists); + const auto cut_width_item = find_label("rpa_cut_width", readinput.input_lists); + ASSERT_NE(cut_mode_item, readinput.input_lists.end()); + ASSERT_NE(cut_radius_item, readinput.input_lists.end()); + ASSERT_NE(cut_width_item, readinput.input_lists.end()); + EXPECT_EQ(TestParameters::input(param).rpa_cut_coulomb, "hard"); + EXPECT_EQ(TestParameters::input(param).rpa_cut_radius, "legacy"); + TestParameters::input(param).rpa_cut_coulomb = "log_erfc"; + TestParameters::input(param).rpa_cut_radius = "inscribed"; + check_item(cut_mode_item->second, param); + check_item(cut_radius_item->second, param); + check_item(cut_width_item->second, param); + TestParameters::input(param).rpa_cut_coulomb = "invalid"; + EXPECT_EXIT(check_item(cut_mode_item->second, param), ::testing::ExitedWithCode(1), ""); + TestParameters::input(param).rpa_cut_radius = "invalid"; + EXPECT_EXIT(check_item(cut_radius_item->second, param), ::testing::ExitedWithCode(1), ""); + TestParameters::input(param).rpa_cut_width = 1.0; + EXPECT_EXIT(check_item(cut_width_item->second, param), ::testing::ExitedWithCode(1), ""); + TestParameters::input(param).rpa_cut_coulomb = "hard"; + TestParameters::input(param).rpa_cut_radius = "legacy"; + TestParameters::input(param).rpa_cut_width = 1.091868692359; + } { // berry_phase auto it = find_label("berry_phase", readinput.input_lists); TestParameters::input(param).berry_phase = true; diff --git a/source/source_io/test_serial/test_read_input.cpp b/source/source_io/test_serial/test_read_input.cpp index 8dd695ec800..801baa2b99d 100644 --- a/source/source_io/test_serial/test_read_input.cpp +++ b/source/source_io/test_serial/test_read_input.cpp @@ -120,6 +120,48 @@ TEST_F(InputDomainReaderTest, PreservesLegacyEmptyValues) EXPECT_TRUE(input.inp.out_pchg.empty()); } +TEST_F(InputDomainReaderTest, RejectsMalformedRpaCutoffTokens) +{ + for (const std::string cut_name : {"rpa_cut_coulomb", "rpa_cut_radius", "rpa_cut_width"}) + { + for (const std::string cut_tokens : {"\nscf_nmax 20\n", "", " invalid\n", " 1 2\n"}) + { + const std::string cut_input = cut_name + cut_tokens; + const std::string cut_output = reject_input(cut_input); + const std::string cut_error = "Invalid INPUT parameter '" + cut_name + "'"; + EXPECT_NE(cut_output.find(cut_error), std::string::npos) << cut_input; + } + } + for (const std::string cut_input : {"rpa_cut_coulomb log_erfc extra\n", + "rpa_cut_radius inscribed extra\n", + "rpa_cut_width 1\n", "rpa_cut_width 1.2501\n", + "rpa_cut_width nan\n", "rpa_cut_width inf\n"}) + { + const std::string cut_output = reject_input(cut_input); + EXPECT_NE(cut_output.find("Invalid INPUT parameter 'rpa_cut_"), std::string::npos) << cut_input; + } +} + +TEST_F(InputDomainReaderTest, AcceptsRpaCutoffChoicesAndWidthUpperBound) +{ + for (const std::string cut_mode : {"hard", "log_erfc"}) + { + for (const std::string cut_radius : {"legacy", "inscribed"}) + { + const std::string cut_input = "rpa_cut_coulomb " + cut_mode + "\nrpa_cut_radius " + + cut_radius + "\nrpa_cut_width 1.25\n"; + write_input(cut_input); + ModuleIO::ReadInput cut_reader(0); + cut_reader.check_ntype_flag = false; + Parameter cut_parameters; + cut_reader.read_parameters(cut_parameters, filename); + EXPECT_EQ(cut_parameters.inp.rpa_cut_coulomb, cut_mode); + EXPECT_EQ(cut_parameters.inp.rpa_cut_radius, cut_radius); + EXPECT_DOUBLE_EQ(cut_parameters.inp.rpa_cut_width, 1.25); + } + } +} + TEST_F(InputDomainReaderTest, RejectsMissingOrdinaryStringValues) { for (const std::string label : {"stru_file", "kpoint_file", "suffix", "dft_functional"}) diff --git a/source/source_lcao/module_ri/CMakeLists.txt b/source/source_lcao/module_ri/CMakeLists.txt index 9f1c903615a..d3ca560b2ec 100644 --- a/source/source_lcao/module_ri/CMakeLists.txt +++ b/source/source_lcao/module_ri/CMakeLists.txt @@ -19,6 +19,7 @@ if(ENABLE_LIBRI) objects lri_cv.cpp conv_coulomb_pot_k.cpp + smooth_cut_coulomb.cpp exx_abfs_ctor_orbs.cpp exx_abfs_io.cpp exx_abfs_jle.cpp diff --git a/source/source_lcao/module_ri/conv_coulomb_pot_k.cpp b/source/source_lcao/module_ri/conv_coulomb_pot_k.cpp index 6903b172154..a0206f47d64 100644 --- a/source/source_lcao/module_ri/conv_coulomb_pot_k.cpp +++ b/source/source_lcao/module_ri/conv_coulomb_pot_k.cpp @@ -1,4 +1,5 @@ #include "conv_coulomb_pot_k.h" +#include "smooth_cut_coulomb.h" #include "source_base/constants.h" #include "source_io/module_parameter/parameter.h" #include "source_basis/module_ao/orb_atomic_lm.h" @@ -162,11 +163,26 @@ namespace Conv_Coulomb_Pot_K case Conv_Coulomb_Pot_K::Coulomb_Type::Fock: { for (const auto& param: param_list.second) { - psik2_ccp = psik2_ccp - + std::stod(param.at("alpha")) - * cal_psi_fock_spencer(orbs.get_psif(), - orbs.get_k_radial(), - std::stod(param.at("Rcut"))); + const auto cut_mode = param.find("rpa_cut_coulomb"); + if (cut_mode != param.end() && cut_mode->second == "log_erfc") + { + const double cut_radius = std::stod(param.at("Rcut")); + const double cut_width = std::stod(param.at("rpa_cut_width")); + const double cut_alpha = std::stod(param.at("alpha")); + for (size_t cut_ik = 0; cut_ik < psik2_ccp.size(); ++cut_ik) + { + psik2_ccp[cut_ik] += cut_alpha * ModuleBase::FOUR_PI * orbs.get_psif()[cut_ik] + * log_erfc_cut_factor(orbs.get_k_radial()[cut_ik], cut_radius, cut_width); + } + } + else + { + psik2_ccp = psik2_ccp + + std::stod(param.at("alpha")) + * cal_psi_fock_spencer(orbs.get_psif(), + orbs.get_k_radial(), + std::stod(param.at("Rcut"))); + } } break; } diff --git a/source/source_lcao/module_ri/exx_lri_detail.cpp b/source/source_lcao/module_ri/exx_lri_detail.cpp index d459fc83790..91dc5582a0f 100644 --- a/source/source_lcao/module_ri/exx_lri_detail.cpp +++ b/source/source_lcao/module_ri/exx_lri_detail.cpp @@ -5,6 +5,9 @@ #include "exx_lri_detail.h" +#include "smooth_cut_coulomb.h" +#include +#include #include "ri_util.h" #include "source_cell/klist.h" @@ -36,12 +39,33 @@ CoulombParam build_center2_cut_coulomb_param(const CoulombParam& coulomb_param, } for (auto& param: param_list.second) { + const auto cut_radius_policy = param.find("rpa_cut_radius"); + if (cut_radius_policy != param.end() && cut_radius_policy->second == "inscribed") + { + const std::array, 3> cut_bvk_lattice = {{ + {{ucell.lat0 * kv.nmp[0] * ucell.a1.x, ucell.lat0 * kv.nmp[0] * ucell.a1.y, ucell.lat0 * kv.nmp[0] * ucell.a1.z}}, + {{ucell.lat0 * kv.nmp[1] * ucell.a2.x, ucell.lat0 * kv.nmp[1] * ucell.a2.y, ucell.lat0 * kv.nmp[1] * ucell.a2.z}}, + {{ucell.lat0 * kv.nmp[2] * ucell.a3.x, ucell.lat0 * kv.nmp[2] * ucell.a3.y, ucell.lat0 * kv.nmp[2] * ucell.a3.z}} + }}; + std::ostringstream cut_radius_text; + cut_radius_text.precision(17); + cut_radius_text << Conv_Coulomb_Pot_K::bvk_inscribed_radius(cut_bvk_lattice); + param["Rcut"] = cut_radius_text.str(); + } auto rcut_it = param.find("Rcut"); if (rcut_it == param.end() || rcut_it->second.empty()) { param["Rcut"] = ModuleBase::GlobalFunc::TO_STRING(fallback_rcut); used_fallback_rcut = true; } + const auto cut_mode = param.find("rpa_cut_coulomb"); + if (cut_mode != param.end()) + { + std::cout << "RPA cut Coulomb: mode=" << cut_mode->second + << " radius=" << param.at("rpa_cut_radius") + << " Rcut_bohr=" << param.at("Rcut") + << " width=" << param.at("rpa_cut_width") << std::endl; + } } } diff --git a/source/source_lcao/module_ri/smooth_cut_coulomb.cpp b/source/source_lcao/module_ri/smooth_cut_coulomb.cpp new file mode 100644 index 00000000000..ea414aecb4c --- /dev/null +++ b/source/source_lcao/module_ri/smooth_cut_coulomb.cpp @@ -0,0 +1,129 @@ +#include "smooth_cut_coulomb.h" + +#include +#include +#include +#include + +namespace Conv_Coulomb_Pot_K +{ +double log_erfc_cut_factor(const double wave_number, + const double cutoff_radius, + const double width_factor) +{ + if (!std::isfinite(wave_number) || wave_number < 0.0 + || !std::isfinite(cutoff_radius) || cutoff_radius <= 0.0 + || !std::isfinite(width_factor) || width_factor <= 1.0 || width_factor > 1.25) + { + throw std::invalid_argument("log-erfc cutoff requires k>=0, Rc>0 and 1 1000000.0) + { + throw std::runtime_error("log-erfc cutoff quadrature exceeds one million panels"); + } + const int cut_panel_count = std::max(128, static_cast(cut_panel_estimate)); + const double cut_step = 16.0 / cut_panel_count; + double cut_integral = 0.0; + for (int cut_panel = 0; cut_panel < cut_panel_count; ++cut_panel) + { + const double cut_coordinate = -8.0 + (cut_panel + 0.5) * cut_step; + const double cut_phase_half = 0.5 * cut_argument * std::exp(cut_log_width * cut_coordinate); + const double cut_sine = std::sin(cut_phase_half); + cut_integral += std::exp(-cut_coordinate * cut_coordinate) * 2.0 * cut_sine * cut_sine; + } + return cut_integral * cut_step / std::sqrt(cut_pi); +} + +double bvk_inscribed_radius(const std::array, 3>& bvk_lattice) +{ + std::array, 3> bvk_cofactor{}; + double bvk_shortest_squared = std::numeric_limits::max(); + for (int bvk_axis = 0; bvk_axis < 3; ++bvk_axis) + { + double bvk_length_squared = 0.0; + for (int cart_axis = 0; cart_axis < 3; ++cart_axis) + { + if (!std::isfinite(bvk_lattice[bvk_axis][cart_axis])) + { + throw std::invalid_argument("non-finite BvK lattice"); + } + bvk_length_squared += bvk_lattice[bvk_axis][cart_axis] * bvk_lattice[bvk_axis][cart_axis]; + const int next_row = (bvk_axis + 1) % 3; + const int last_row = (bvk_axis + 2) % 3; + const int next_col = (cart_axis + 1) % 3; + const int last_col = (cart_axis + 2) % 3; + bvk_cofactor[bvk_axis][cart_axis] + = bvk_lattice[next_row][next_col] * bvk_lattice[last_row][last_col] + - bvk_lattice[next_row][last_col] * bvk_lattice[last_row][next_col]; + } + bvk_shortest_squared = std::min(bvk_shortest_squared, bvk_length_squared); + } + double bvk_determinant = 0.0; + for (int cart_axis = 0; cart_axis < 3; ++cart_axis) + { + bvk_determinant += bvk_lattice[0][cart_axis] * bvk_cofactor[0][cart_axis]; + } + if (bvk_determinant == 0.0 || !std::isfinite(bvk_determinant) || bvk_shortest_squared <= 0.0) + { + throw std::invalid_argument("singular BvK lattice"); + } + // For d=n*B, |n_i| <= |d|*||(B^-1)[:,i]||. These bounds include + // every vector shorter than the shortest supplied basis vector. + std::array bvk_bounds{}; + double bvk_candidate_count = 1.0; + for (int bvk_axis = 0; bvk_axis < 3; ++bvk_axis) + { + double bvk_inverse_norm_squared = 0.0; + for (int cart_axis = 0; cart_axis < 3; ++cart_axis) + { + const double bvk_inverse_entry = bvk_cofactor[bvk_axis][cart_axis] / bvk_determinant; + bvk_inverse_norm_squared += bvk_inverse_entry * bvk_inverse_entry; + } + const double bvk_bound = std::floor(std::sqrt(bvk_shortest_squared * bvk_inverse_norm_squared) + 1.e-10); + if (!std::isfinite(bvk_bound) || bvk_bound > 1000000.0) + { + throw std::runtime_error("BvK shortest-vector search needs a reduced lattice"); + } + bvk_bounds[bvk_axis] = static_cast(bvk_bound); + bvk_candidate_count *= 2.0 * bvk_bounds[bvk_axis] + 1.0; + } + if (bvk_candidate_count > 10000000.0) + { + throw std::runtime_error("BvK shortest-vector search exceeds ten million candidates"); + } + for (int first_cell = -bvk_bounds[0]; first_cell <= bvk_bounds[0]; ++first_cell) + { + for (int second_cell = -bvk_bounds[1]; second_cell <= bvk_bounds[1]; ++second_cell) + { + for (int third_cell = -bvk_bounds[2]; third_cell <= bvk_bounds[2]; ++third_cell) + { + if (first_cell == 0 && second_cell == 0 && third_cell == 0) { continue; } + double bvk_distance_squared = 0.0; + for (int cart_axis = 0; cart_axis < 3; ++cart_axis) + { + const double bvk_displacement = first_cell * bvk_lattice[0][cart_axis] + + second_cell * bvk_lattice[1][cart_axis] + third_cell * bvk_lattice[2][cart_axis]; + bvk_distance_squared += bvk_displacement * bvk_displacement; + } + bvk_shortest_squared = std::min(bvk_shortest_squared, bvk_distance_squared); + } + } + } + return 0.5 * std::sqrt(bvk_shortest_squared); +} +} diff --git a/source/source_lcao/module_ri/smooth_cut_coulomb.h b/source/source_lcao/module_ri/smooth_cut_coulomb.h new file mode 100644 index 00000000000..4b070872b23 --- /dev/null +++ b/source/source_lcao/module_ri/smooth_cut_coulomb.h @@ -0,0 +1,15 @@ +#ifndef ABACUS_SMOOTH_CUT_COULOMB_H +#define ABACUS_SMOOTH_CUT_COULOMB_H + +#include + +namespace Conv_Coulomb_Pot_K +{ +// k^2 v(k)/(4*pi) for f(r)/r, f=erfc(log(r/Rc)/log(width))/2. +double log_erfc_cut_factor(double wave_number, double cutoff_radius, double width_factor); + +// Half the shortest nonzero superlattice vector; row vectors are in Bohr. +double bvk_inscribed_radius(const std::array, 3>& bvk_lattice); +} + +#endif diff --git a/source/source_lcao/module_ri/unittests/CMakeLists.txt b/source/source_lcao/module_ri/unittests/CMakeLists.txt index 887eaa63c1b..17c4693fd20 100644 --- a/source/source_lcao/module_ri/unittests/CMakeLists.txt +++ b/source/source_lcao/module_ri/unittests/CMakeLists.txt @@ -1,6 +1,9 @@ abacus_disable_feature_definitions(__MLALGO) abacus_disable_feature_definitions(__CUDA) abacus_disable_feature_definitions(__ROCM) +addtest( + TARGET MODULE_RI_smooth_cut_coulomb + SOURCES test_smooth_cut_coulomb.cpp ../smooth_cut_coulomb.cpp) if(ENABLE_LCAO) addtest( TARGET diff --git a/source/source_lcao/module_ri/unittests/test_smooth_cut_coulomb.cpp b/source/source_lcao/module_ri/unittests/test_smooth_cut_coulomb.cpp new file mode 100644 index 00000000000..f14a803d747 --- /dev/null +++ b/source/source_lcao/module_ri/unittests/test_smooth_cut_coulomb.cpp @@ -0,0 +1,34 @@ +#include "../smooth_cut_coulomb.h" +#include +#include +#include + +TEST(SmoothCutCoulomb, GammaMomentAndSharpLimit) +{ + const double cut_width = 1.091868692359; + const double cut_sigma = std::log(cut_width); + const double cut_z = 1.e-7; + EXPECT_NEAR(Conv_Coulomb_Pot_K::log_erfc_cut_factor(cut_z, 1.0, cut_width) / (cut_z * cut_z), + 0.5 * std::exp(cut_sigma * cut_sigma), 1.e-13); + EXPECT_DOUBLE_EQ(Conv_Coulomb_Pot_K::log_erfc_cut_factor(0.0, 32.0, cut_width), 0.0); + EXPECT_NEAR(Conv_Coulomb_Pot_K::log_erfc_cut_factor(0.2, 32.0, 1.00000001), + 1.0 - std::cos(6.4), 1.e-13); + for (int cut_index = 1; cut_index <= 100; ++cut_index) + { + const double cut_factor = Conv_Coulomb_Pot_K::log_erfc_cut_factor(cut_index, 20.0, cut_width); + EXPECT_GE(cut_factor, 0.0); + EXPECT_LE(cut_factor, 2.0 + 1.e-14); + } +} + +TEST(SmoothCutCoulomb, InscribedRadiusUsesShortestCombination) +{ + const std::array, 3> cut_cubic = {{{{8., 0., 0.}}, {{0., 10., 0.}}, {{0., 0., 12.}}}}; + const std::array, 3> cut_fcc = {{{{0., 4., 4.}}, {{4., 0., 4.}}, {{4., 4., 0.}}}}; + const std::array, 3> cut_skew = {{{{10., 0., 0.}}, {{9.9, 0.1, 0.}}, {{0., 0., 5.}}}}; + EXPECT_DOUBLE_EQ(Conv_Coulomb_Pot_K::bvk_inscribed_radius(cut_cubic), 4.0); + EXPECT_NEAR(Conv_Coulomb_Pot_K::bvk_inscribed_radius(cut_fcc), std::sqrt(8.0), 1.e-14); + EXPECT_NEAR(Conv_Coulomb_Pot_K::bvk_inscribed_radius(cut_skew), std::sqrt(0.02) / 2., 1.e-14); + EXPECT_THROW(Conv_Coulomb_Pot_K::log_erfc_cut_factor(1., 1., 1.), std::invalid_argument); + EXPECT_THROW(Conv_Coulomb_Pot_K::log_erfc_cut_factor(-1., 1., 1.1), std::invalid_argument); +}