mohanchen · Oct 13, 2026 · Sep 30, 2026 · Sep 30, 2026 · Sep 30, 2026 · Sep 30, 2026
diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
@@ -206,6 +206,7 @@
     - [rpa](#rpa)
     - [rpa\_out\_vel](#rpa_out_vel)
     - [rpa\_outdir](#rpa_outdir)
+    - [out\_librpa\_ver](#out_librpa_ver)
     - [out\_pchg](#out_pchg)
     - [out\_wfc\_norm](#out_wfc_norm)
     - [out\_wfc\_re\_im](#out_wfc_re_im)
@@ -2327,7 +2328,7 @@
 
 - **Type**: Boolean
 - **Availability**: *[`basis_type`](#basis_type) in [lcao, lcao_in_pw]*
-- **Description**: Whether to print the upper triangular part of the exchange-correlation matrices in Kohn-Sham orbital representation: for each k point into files in the directory OUT.i_nao.txt, where {suffix}/vxc_out.dat. If EXX is calculated, the local and EXX part of band energy will also be printed in OUT.{suffix}/vxc_exx_out.dat, respectively. All the vxc_out.dat files contains 3 integers (nk, nspin, nband) followed by nk*nspin*nband lines of energy Hartree and eV.
+- **Description**: Write exchange-correlation matrices in the Kohn-Sham orbital basis and their diagonal expectation values. The diagonal energies are saved in OUT.${suffix}/vxc.txt. When EXX is calculated, its local and EXX contributions are also saved as vxc_local.txt and vxc_exx.txt in the same directory. Each energy file begins with nk, nspin, and nbands, followed by nk*nspin*nbands rows of energies in Hartree and eV.
 
   > Note: In the 3.10-LTS version, the file name is k-$k-Vxc and so on.
 - **Default**: False
@@ -2456,6 +2457,13 @@
 - **Description**: The directory to save files for LibRPA.
 - **Default**: "OUT.librpa"
 
+### out_librpa_ver
+
+- **Type**: Integer
+- **Availability**: *[`basis_type`](#basis_type)==lcao*
+- **Description**: Select the LibRPA file protocol: 0 for legacy text output and 1 for reader-v1 output. Reader-v1 supports symmetry-reduced Monkhorst-Pack grids and requires a uniform Monkhorst-Pack KPT grid.
+- **Default**: 0
+
 ### out_pchg
 
 - **Type**: String

diff --git a/docs/parameters.yaml b/docs/parameters.yaml
@@ -3734,7 +3734,7 @@ parameters:
     category: Output information
     type: Boolean
     description: |
-      Whether to print the upper triangular part of the exchange-correlation matrices in Kohn-Sham orbital representation: for each k point into files in the directory OUT.i_nao.txt, where {suffix}/vxc_out.dat. If EXX is calculated, the local and EXX part of band energy will also be printed in OUT.{suffix}/vxc_exx_out.dat, respectively. All the vxc_out.dat files contains 3 integers (nk, nspin, nband) followed by nk*nspin*nband lines of energy Hartree and eV.
+      Write exchange-correlation matrices in the Kohn-Sham orbital basis and their diagonal expectation values. The diagonal energies are saved in OUT.${suffix}/vxc.txt. When EXX is calculated, its local and EXX contributions are also saved as vxc_local.txt and vxc_exx.txt in the same directory. Each energy file begins with nk, nspin, and nbands, followed by nk*nspin*nbands rows of energies in Hartree and eV.
 
       [NOTE] In the 3.10-LTS version, the file name is k-$k-Vxc and so on.
     default_value: "False"
@@ -3876,6 +3876,14 @@ parameters:
     default_value: "\"OUT.librpa\""
     unit: ""
     availability: basis_type==lcao
+  - name: out_librpa_ver
+    category: Output information
+    type: Integer
+    description: |
+      Select the LibRPA file protocol: 0 for legacy text output and 1 for reader-v1 output. Reader-v1 supports symmetry-reduced Monkhorst-Pack grids and requires a uniform Monkhorst-Pack KPT grid.
+    default_value: "0"
+    unit: ""
+    availability: basis_type==lcao
   - name: out_pchg
     category: Output information
     type: String

diff --git a/source/Makefile.Objects b/source/Makefile.Objects
@@ -858,6 +858,24 @@ OBJS_MODULE_RI=conv_coulomb_pot_k.o\
       gaussian_abfs.o\
       exx_lri_detail.o\
       ri_util.o\
+      rpa_lri.o\
+      rpa_lri_detail.o\
+      rpa_lri_band_io.o\
+      rpa_lri_mpi_layout.o\
+      rpa_lri_wfc_mpi.o\
+      rpa_lri_bz_io.o\
+      rpa_lri_coulomb.o\
+      rpa_lri_coulomb_io.o\
+      rpa_lri_overlap.o\
+      rpa_lri_overlap_io.o\
+      rpa_lri_output.o\
+      rpa_lri_bz.o\
+      rpa_lri_sym.o\
+      rpa_lri_coeff_io.o\
+      rpa_lri_velocity.o\
+      rpa_lri_basis.o\
+      rpa_lri_eigenvector_v1.o\
+      rpa_lri_eigenvector.o\
       singular_value.o\
 
 OBJS_BSE=esolver_lr_lcao_bse.o\

diff --git a/source/source_cell/klist.cpp b/source/source_cell/klist.cpp
@@ -87,6 +87,13 @@ void K_Vectors::set(const UnitCell& ucell,
                             this->kvec_c_full);
 
 
+    // initialize ibz_index
+    this->ibz_index.resize(this->nkstot_nospin);
+    for (int ik = 0; ik < this->nkstot_nospin; ik++)
+    {
+        this->ibz_index[ik] = ik;
+    }
+
     // (2)
     // reduce kpoints to IBZ according to symmetry operations
     if (use_ibz)
@@ -140,13 +147,6 @@ void K_Vectors::set(const UnitCell& ucell,
     // set the k vectors for the up and down spin
     this->set_kup_and_kdw(ofs);
 
-    // initialize ibz_index
-    this->ibz_index.resize(this->nkstot_nospin);
-    for (int ik = 0; ik < this->nkstot_nospin; ik++)
-    {
-        this->ibz_index[ik] = ik;
-    }
-
     // get ik2iktot: map local k indices to global indices in the pool
     KListIO::build_ik2iktot(this->para_k.my_pool,
                             this->para_k.startk_pool,
@@ -203,6 +203,24 @@ void K_Vectors::renew(const int& kpoint_number)
     return;
 }
 
+bool K_Vectors::is_uniform_mp() const
+{
+    if (!this->is_mp || this->nkstot_nospin <= 0)
+    {
+        return false;
+    }
+    long long full_kpoints = 1;
+    for (const int size : this->nmp)
+    {
+        if (size <= 0 || full_kpoints > this->nkstot_nospin / size)
+        {
+            return false;
+        }
+        full_kpoints *= size;
+    }
+    return full_kpoints == this->nkstot_nospin;
+}
+
 // Read the KPT file, which contains K-point coordinates, weights, and grid size information
 // Generate K-point grid according to different parameters of the KPT file
 bool K_Vectors::read_kpoints(const UnitCell& ucell,
@@ -663,6 +681,9 @@ void K_Vectors::mpi_k(std::ofstream& ofs_running, const int my_rank, const int m
 
     Parallel_Common::bcast_int(this->nmp, 3);
 
+    // The MP flag is set while rank 0 reads KPT and is consumed on every rank.
+    Parallel_Common::bcast_bool(this->is_mp);
+
     this->kl_segids.resize(this->nkstot);
     Parallel_Common::bcast_int(this->kl_segids.data(), this->nkstot);
 

diff --git a/source/source_cell/klist.h b/source/source_cell/klist.h
@@ -132,6 +132,8 @@ class K_Vectors : public ModuleCell::ReciprocalGrid
         return is_mp;
     }
 
+    bool is_uniform_mp() const;
+
     std::vector<int> ik2iktot; ///<[nks] map ik to the global index of k points
     std::vector<int> ibz_index; ///< map k points (before symmetry reduction) to irreducible k-points
 
@@ -457,4 +459,4 @@ class K_Vectors : public ModuleCell::ReciprocalGrid
     void mpi_k(std::ofstream& ofs_running, const int my_rank, const int my_pool);
 #endif
 };
-#endif // KVECT_H
+#endif // KVECT_H
diff --git a/source/source_cell/klist_io.cpp b/source/source_cell/klist_io.cpp
@@ -271,6 +271,11 @@ void pack_kpts(const std::vector<int>& isk,
         kvec_d_aux[3 * ik] = kvec_d[ik].x;
         kvec_d_aux[3 * ik + 1] = kvec_d[ik].y;
         kvec_d_aux[3 * ik + 2] = kvec_d[ik].z;
+    }
+    // The unreduced mesh is global and can be larger than the SCF list.
+    const std::size_t full_count = std::min(kvec_c_full.size(), kvec_c_full_aux.size() / 3);
+    for (std::size_t ik = 0; ik < full_count; ++ik)
+    {
         kvec_c_full_aux[3 * ik] = kvec_c_full[ik].x;
         kvec_c_full_aux[3 * ik + 1] = kvec_c_full[ik].y;
         kvec_c_full_aux[3 * ik + 2] = kvec_c_full[ik].z;
@@ -300,12 +305,16 @@ void unpack_kpts(const std::vector<int>& isk_aux,
         kvec_d[i].x = kvec_d_aux[k_index * 3];
         kvec_d[i].y = kvec_d_aux[k_index * 3 + 1];
         kvec_d[i].z = kvec_d_aux[k_index * 3 + 2];
-        kvec_c_full[i].x = kvec_c_full_aux[k_index * 3];
-        kvec_c_full[i].y = kvec_c_full_aux[k_index * 3 + 1];
-        kvec_c_full[i].z = kvec_c_full_aux[k_index * 3 + 2];
         wk[i] = wk_aux[k_index];
         isk[i] = isk_aux[k_index];
     }
+    kvec_c_full.resize(kvec_c_full_aux.size() / 3);
+    for (std::size_t ik = 0; ik < kvec_c_full.size(); ++ik)
+    {
+        kvec_c_full[ik].x = kvec_c_full_aux[3 * ik];
+        kvec_c_full[ik].y = kvec_c_full_aux[3 * ik + 1];
+        kvec_c_full[ik].z = kvec_c_full_aux[3 * ik + 2];
+    }
 }
 
 void bcast_kstars(std::vector<std::map<int, ModuleBase::Vector3<double>>>& kstars,

diff --git a/source/source_cell/test/CMakeLists.txt b/source/source_cell/test/CMakeLists.txt
@@ -95,7 +95,7 @@ AddTest(
 AddTest(
   TARGET MODULE_CELL_klist_test
   LIBS base device symmetry
-  SOURCES klist_test.cpp ../klist.cpp ../klist_io.cpp ../parallel_kpoints.cpp ../reciprocal_grid.cpp
+  SOURCES klist_test.cpp test_klist_io.cpp ../klist.cpp ../klist_io.cpp ../parallel_kpoints.cpp ../reciprocal_grid.cpp
 )
 
 AddTest(

diff --git a/source/source_cell/test/klist_test.cpp b/source/source_cell/test/klist_test.cpp
@@ -881,6 +881,7 @@ TEST_F(KlistTest, IbzKpointIsMP)
     kv->read_kpoints_for_testing(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank);
     EXPECT_EQ(kv->get_nkstot(), 512);
     EXPECT_TRUE(kv->is_mp);
+    EXPECT_TRUE(kv->is_uniform_mp());
     // calculate ibz_kpoint
     std::string skpt;
     ModuleSymmetry::Symmetry::symm_flag = 0;
@@ -915,6 +916,7 @@ TEST_F(KlistTest, IbzKpointCustomWeights)
         kv_test1.read_kpoints_for_testing(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank);
         EXPECT_EQ(kv_test1.get_nkstot(), 5);
         EXPECT_FALSE(kv_test1.is_mp); // Should be non-MP
+        EXPECT_FALSE(kv_test1.is_uniform_mp());
 
         // Store original weights before IBZ reduction
         std::vector<double> original_weights = kv_test1.wk;
@@ -944,6 +946,7 @@ TEST_F(KlistTest, IbzKpointCustomWeights)
         kv_test2.read_kpoints_for_testing(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank);
         EXPECT_EQ(kv_test2.get_nkstot(), 5);
         EXPECT_FALSE(kv_test2.is_mp); // Should be non-MP
+        EXPECT_FALSE(kv_test2.is_uniform_mp());
 
         // Verify custom weights were read correctly
         EXPECT_DOUBLE_EQ(kv_test2.wk[0], 0.1);
@@ -1001,6 +1004,7 @@ TEST_F(KlistTest, IbzKpointCustomWeights)
         kv_test3.read_kpoints_for_testing(ucell, k_file, gamma_only_local, kspacing, kmesh_type, koffset, ofs_running, ofs_warning, my_rank);
         EXPECT_EQ(kv_test3.get_nkstot(), 512);
         EXPECT_TRUE(kv_test3.is_mp); // Should be MP
+        EXPECT_TRUE(kv_test3.is_uniform_mp());
 
         // Apply IBZ reduction
         std::string skpt;

diff --git a/source/source_cell/test/klist_test_para.cpp b/source/source_cell/test/klist_test_para.cpp
@@ -208,6 +208,7 @@ TEST_F(KlistParaTest, Set)
     const std::string kmesh_type = "gamma";
     const double koffset[3] = {0.0, 0.0, 0.0};
     kv->set(ucell, symm, k_file, /*nspin_in*/ 1, ucell.G, ucell.latvec, GlobalV::ofs_running, GlobalV::ofs_warning, use_ibz, global_out_dir, gamma_only_local, kspacing, kmesh_type, koffset);
+    EXPECT_TRUE(kv->get_is_mp());
     EXPECT_EQ(kv->get_nkstot(), 35);
     EXPECT_EQ(kv->get_nkstot_nospin(), 512);
     EXPECT_GT(kv->get_nkstot_nospin(), kv->get_nkstot());

diff --git a/source/source_cell/test/test_klist_io.cpp b/source/source_cell/test/test_klist_io.cpp
@@ -0,0 +1,39 @@
+#include "../klist_io.h"
+#include <gtest/gtest.h>
+
+TEST(KListIO, PreserveFullMeshAcrossPools)
+{
+    const std::vector<int> spins = {0, 0};
+    const std::vector<double> weights = {1.0 / 3, 2.0 / 3};
+    const std::vector<ModuleBase::Vector3<double>> irreducible = {{0, 0, 0}, {0, 0, 1.0 / 3}};
+    const std::vector<ModuleBase::Vector3<double>> full = {{0, 0, 0}, {0, 0, 1.0 / 3}, {0, 0, 2.0 / 3}};
+    auto full_with_spin_capacity = full;
+    full_with_spin_capacity.resize(6);
+    std::vector<int> spin_buffer(2);
+    std::vector<double> weight_buffer(2);
+    std::vector<double> cart_buffer(6);
+    std::vector<double> frac_buffer(6);
+    std::vector<double> full_buffer(9, -9.0);
+    KListIO::pack_kpts(spins, weights, irreducible, irreducible, full_with_spin_capacity, 2,
+                       spin_buffer, weight_buffer, cart_buffer, frac_buffer, full_buffer);
+    EXPECT_DOUBLE_EQ(full_buffer[8], 2.0 / 3);
+
+    std::vector<int> local_spin(1);
+    std::vector<double> local_weight(1);
+    std::vector<ModuleBase::Vector3<double>> local_cart(1);
+    std::vector<ModuleBase::Vector3<double>> local_frac(1);
+    std::vector<ModuleBase::Vector3<double>> local_full(1);
+    // A second pool holds only the second SCF representative, but must retain
+    // every original full-grid point in its original order.
+    KListIO::unpack_kpts(spin_buffer, weight_buffer, cart_buffer, frac_buffer, full_buffer,
+                         1, 1, local_spin, local_weight, local_cart, local_frac, local_full);
+    EXPECT_DOUBLE_EQ(local_cart[0].z, 1.0 / 3);
+    EXPECT_DOUBLE_EQ(local_weight[0], 2.0 / 3);
+    ASSERT_EQ(local_full.size(), full.size());
+    for (std::size_t ik = 0; ik < full.size(); ++ik)
+    {
+        EXPECT_DOUBLE_EQ(local_full[ik].x, full[ik].x);
+        EXPECT_DOUBLE_EQ(local_full[ik].y, full[ik].y);
+        EXPECT_DOUBLE_EQ(local_full[ik].z, full[ik].z);
+    }
+}
diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp
@@ -22,7 +22,6 @@
 #include "source_base/module_parallel/para_world.h"
 #include "source_base/module_parallel/para_tag.h"
 #include "source_base/module_parallel/para_bridge.h"
-
 namespace ModuleESolver
 {
 

diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp
@@ -1,6 +1,8 @@
 #include "esolver_ks_lcao.h"
 #include "source_basis/module_pw/pw_basis_big.h" // use PW_Basis_Big
+#include "source_base/tool_quit.h"
 #include "source_base/module_external/blacs_connector.h"
+#include "source_cell/klist.h"
 #include "source_cell/module_neighbor/sltk_atom_arrange.h"
 #include "source_estate/elecstate_tools.h"
 #include "source_lcao/module_deltaspin/spin_constrain.h"
@@ -24,6 +26,7 @@
 #include "source_estate/module_charge/chg_routine.h" // use charge mixing, mohan add 20251006
 #include "source_estate/module_dm/dm_routine.h" // init dm from electronic wave functions
 #include "source_io/module_restart/restart.h" // GlobalC::restart for load_exx_flag
+#include "source_io/module_parameter/input_parameter.h"
 #include "source_io/module_ctrl/ctrl_runner_lcao.h" // use ctrl_runner_lcao()
 #include "source_io/module_ctrl/ctrl_iter_lcao.h" // use ctrl_iter_lcao()
 #include "source_io/module_ctrl/ctrl_scf_lcao.h" // use ctrl_scf_lcao()
@@ -36,6 +39,22 @@
 namespace ModuleESolver
 {
 
+#ifdef __EXX
+namespace
+{
+void validate_librpa_v1_kpoints(const Input_para& input, const K_Vectors& kv)
+{
+    // Reader-v1 writes the complete Monkhorst-Pack mesh for LibRPA.
+    if (input.rpa && input.out_librpa_ver == 1 && !kv.is_uniform_mp())
+    {
+        ModuleBase::WARNING_QUIT(
+            "ESolver_KS_LCAO::before_all_runners",
+            "out_librpa_ver=1 requires a uniform Monkhorst-Pack k-point grid.");
+    }
+}
+} // namespace
+#endif
+
 template <typename TK, typename TR>
 ESolver_KS_LCAO<TK, TR>::ESolver_KS_LCAO()
 {
@@ -65,6 +84,10 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(BaseCell& basecell, const Input
     // 1) before_all_runners in ESolver_KS (includes init_general_exx_info)
     ESolver_KS::before_all_runners(ucell, inp);
 
+#ifdef __EXX
+    validate_librpa_v1_kpoints(inp, this->kv);
+#endif
+
     // 2) init full Exx_Info for LCAO (includes info_ri, info_opt_abfs, info_lip)
     init_exx_info(this->exx_info_, inp);
 

diff --git a/source/source_io/module_ctrl/ctrl_scf_lcao.cpp b/source/source_io/module_ctrl/ctrl_scf_lcao.cpp
@@ -683,7 +683,14 @@ void ModuleIO::ctrl_scf_lcao(UnitCell& ucell,
     //------------------------------------------------------------------
     if (inp.rpa)
     {
-        RPA_LRI<TK, double> rpa_lri_double(exx_info.info_ri);
+        const ModuleRI::RpaLriRuntime rpa_runtime = {
+            inp,
+            global_out_dir,
+            pv.get_global_row_size(),
+            GlobalV::MY_RANK,
+            GlobalV::NPROC,
+            GlobalV::ofs_running};
+        RPA_LRI<TK, double> rpa_lri_double(exx_info.info_ri, rpa_runtime);
         rpa_lri_double.postSCF(ucell, MPI_COMM_WORLD, *dm, pelec, kv, orb, pv, *psi);
         if (inp.rpa_out_vel)
             rpa_lri_double.out_velocity(ucell, gd, two_center_bundle, pv, *psi, pelec);

diff --git a/source/source_io/module_hs/vxc_op_lip.cpp b/source/source_io/module_hs/vxc_op_lip.cpp
@@ -219,8 +219,8 @@ void write_Vxc_LIP(int nspin,
             assert(e_orb.size() == kv.get_nks());
             const int nk = kv.get_nks() / nspin0;
             std::ofstream ofs;
-            const std::string out_name = (label == "") ? "out.dat" : label + "_out.dat";
-            ofs.open(global_out_dir + "vxc_" + out_name,
+            const std::string out_name = (label == "") ? ".txt" : "_" + label + ".txt";
+            ofs.open(global_out_dir + "vxc" + out_name,
                 app ? std::ios::app : std::ios::out);
             ofs << nk << "\n" << nspin0 << "\n" << nbands << "\n";
             ofs << std::scientific << std::setprecision(16);

diff --git a/source/source_io/module_hs/vxc_op_tools.cpp b/source/source_io/module_hs/vxc_op_tools.cpp
@@ -110,8 +110,10 @@ void write_orb_energy(const K_Vectors& kv,
     assert(e_orb.size() == kv.get_nks());
     const int nk = kv.get_nks() / nspin0;
     std::ofstream ofs;
-    const std::string out_name = (label == "") ? "out.dat" : label + "_out.dat";
-    ofs.open(global_out_dir + term + "_" + out_name,
+    const std::string out_name = term == "vxc"
+        ? (label.empty() ? "vxc.txt" : "vxc_" + label + ".txt")
+        : term + "_" + (label.empty() ? "out.dat" : label + "_out.dat");
+    ofs.open(global_out_dir + out_name,
         app ? std::ios::app : std::ios::out);
     ofs << nk << "\n" << nspin0 << "\n" << nbands << "\n";
     ofs << std::scientific << std::setprecision(16);

diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h
@@ -648,6 +648,7 @@ struct Input_para
     double shrink_LU_inv_thr = 1e-6; ///< threshold to get inverse of overlap matrix by LU decomposition in auxiliary basis representation
     bool out_ri_cv = false;          ///< Whether to output the coefficient tensor C and ABFs-representation Coulomb matrix V
     bool out_unshrinked_v = false;   ///< whether to output the large Vq matrix in unshrinked auxiliary basis
+    int out_librpa_ver = 0; ///< 0: legacy ABACUS LibRPA text files; 1: LibRPA reader-v1 binary files
     bool exx_coul_moment = false;    ///< whether to use moment method for Coulomb calculation
     bool exx_rotate_abfs = false;    ///< whether to rotate auxiliary basis for Coulomb calculation
     double exx_multip_moments_threshold = 1e-10; ///< threshold to screen multipole moments in Coulomb calculation