Add RPA smooth coulomb cutoff for V_cut

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\

source/source_hamilt/module_xc/exx_info.cpp

@@ -8,6 +8,8 @@
 #include <algorithm>
 #include <cassert>
 #include <string>
+#include <sstream>
+#include <stdexcept>
 #include <vector>
 
 //----------------------------------------------------------
@@ -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;

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

source/source_hamilt/module_xc/test/test_exx_info.cpp

@@ -0,0 +1,41 @@
+#include "../exx_info.h"
+#include "source_io/module_parameter/input_parameter.h"
+#include <gtest/gtest.h>
+#include <stdexcept>
+
+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);
+
+    cut_input.dft_functional = "pbe";
+    cut_input.rpa_cut_coulomb = "log_erfc";
+    Exx_Info cut_pbe_info;
+    ASSERT_NO_THROW(init_exx_info(cut_pbe_info, cut_input));
+    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");
+}

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

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";

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;

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"})

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

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;
             }

source/source_lcao/module_ri/exx_lri_detail.cpp

@@ -5,6 +5,9 @@
 
 
 #include "exx_lri_detail.h"
+#include "smooth_cut_coulomb.h"
+#include <iostream>
+#include <sstream>
 
 #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<std::array<double, 3>, 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;
+            }
         }
     }