diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt
index 9f7ab0196cc..f613fbe3e01 100644
--- a/source/CMakeLists.txt
+++ b/source/CMakeLists.txt
@@ -572,6 +572,7 @@ add_subdirectory(source_esolver)
 add_subdirectory(source_io)
 add_subdirectory(source_relax)
 add_subdirectory(source_lcao/module_ri)
+add_subdirectory(source_lcao/module_rpa)
 add_subdirectory(source_io/module_parameter)
 add_subdirectory(source_lcao/module_lr)
 add_subdirectory(source_lcao/module_rdmft) # add by jghan
@@ -800,6 +801,9 @@ if(ENABLE_MLALGO)
 endif()
 if(ENABLE_LIBRI)
   target_link_libraries(${ABACUS_BIN_NAME} PRIVATE ri module_exx_symmetry)
+  if(ENABLE_LCAO)
+    target_link_libraries(${ABACUS_BIN_NAME} PRIVATE rpa librpa_io)
+  endif()
 endif()
 if(ENABLE_JSON)
   target_link_libraries(${ABACUS_BIN_NAME} PRIVATE json_output)
diff --git a/source/Makefile b/source/Makefile
index 6818478bbae..afd5d9364eb 100644
--- a/source/Makefile
+++ b/source/Makefile
@@ -144,6 +144,9 @@ endif
 ifdef LIBRI_DIR
     INCLUDES += -I${LIBRI_DIR}/include
     INCLUDES += -I${LIBCOMM_DIR}/include
+    INCLUDES += -I${ABACUS_ROOT}source_lcao/module_ri
+    INCLUDES += -I${ABACUS_ROOT}source_lcao/module_rpa
+    INCLUDES += -I${ABACUS_ROOT}source_lcao/module_rpa/module_librpa_io
     HONG                += -DUSE_LIBRI
     HONG += -D__EXX -DEXX_H_COMM=2 -DEXX_DM=3 -DTEST_EXX_LCAO=0 -DTEST_EXX_RADIAL=1
 endif
@@ -241,7 +244,7 @@ HONG += -D__SELINV -DMETIS
 
 ifeq ($(MPI), ON)
     ifdef LIBRI_DIR
-        OBJS_ABACUS += $(OBJS_MODULE_RI) $(OBJS_BSE) ${OBJS_SRC_RI}
+        OBJS_ABACUS += $(OBJS_MODULE_RI) $(OBJS_MODULE_RPA) $(OBJS_BSE) ${OBJS_SRC_RI}
     endif
 endif
 
diff --git a/source/Makefile.Objects b/source/Makefile.Objects
index 34b39f96da5..4956896a05a 100644
--- a/source/Makefile.Objects
+++ b/source/Makefile.Objects
@@ -76,6 +76,8 @@ VPATH=./src_global:\
 ./src_ri:\
 ./source_lcao/module_ri:\
 ./source_lcao/module_ri/module_exx_symmetry:\
+./source_lcao/module_rpa:\
+./source_lcao/module_rpa/module_librpa_io:\
 ./source_lcao/module_bse:\
 ./source_io/module_parameter:\
 ./source_io/module_efield:\
@@ -857,9 +859,12 @@ OBJS_MODULE_RI=conv_coulomb_pot_k.o\
       gaussian_abfs.o\
       exx_lri_detail.o\
       ri_util.o\
-      rpa_lri.o\
+      singular_value.o\
+
+OBJS_MODULE_RPA=rpa_lri.o\
       rpa_lri_detail.o\
       rpa_lri_band_io.o\
+      rpa_lri_coeff_v1.o\
       rpa_lri_mpi_layout.o\
       rpa_lri_wfc_mpi.o\
       rpa_lri_bz_io.o\
@@ -870,12 +875,14 @@ OBJS_MODULE_RI=conv_coulomb_pot_k.o\
       rpa_lri_output.o\
       rpa_lri_bz.o\
       rpa_lri_sym.o\
+      librpa_bz_sampling.o\
+      librpa_stru_symmetry.o\
+      librpa_stru_units.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\
     hamilt_bse.o\
diff --git a/source/source_io/module_ctrl/ctrl_scf_lcao.cpp b/source/source_io/module_ctrl/ctrl_scf_lcao.cpp
index b7533f0325d..6c21487e0a4 100644
--- a/source/source_io/module_ctrl/ctrl_scf_lcao.cpp
+++ b/source/source_io/module_ctrl/ctrl_scf_lcao.cpp
@@ -37,7 +37,7 @@
 #endif
 #ifdef __EXX
 #include "source_lcao/module_ri/exx_lri_interface.h" // use EXX codes
-#include "source_lcao/module_ri/rpa_lri.h"           // use RPA code
+#include "source_lcao/module_rpa/rpa_lri.h"          // use RPA code
 #endif
 #include "../module_qo/to_qo.h"                // use toQO
 #include "source_lcao/module_rdmft/rdmft.h" // use RDMFT codes
diff --git a/source/source_lcao/module_ri/CMakeLists.txt b/source/source_lcao/module_ri/CMakeLists.txt
index 8b0e6e39494..aa4bdb6059c 100644
--- a/source/source_lcao/module_ri/CMakeLists.txt
+++ b/source/source_lcao/module_ri/CMakeLists.txt
@@ -23,25 +23,6 @@ if (ENABLE_LIBRI)
         singular_value.cpp
         exx_lri_detail.cpp
         ri_util.cpp
-        rpa_lri.cpp
-        rpa_lri_detail.cpp
-        rpa_lri_band_io.cpp
-        rpa_lri_coeff_v1.cpp
-        rpa_lri_bz_io.cpp
-        rpa_lri_mpi_layout.cpp
-        rpa_lri_wfc_mpi.cpp
-        rpa_lri_coulomb.cpp
-        rpa_lri_coulomb_io.cpp
-        rpa_lri_overlap.cpp
-        rpa_lri_overlap_io.cpp
-        rpa_lri_output.cpp
-        rpa_lri_bz.cpp
-        rpa_lri_sym.cpp
-        rpa_lri_coeff_io.cpp
-        rpa_lri_velocity.cpp
-        rpa_lri_basis.cpp
-        rpa_lri_eigenvector_v1.cpp
-        rpa_lri_eigenvector.cpp
     )
   endif()
   add_library(
diff --git a/source/source_lcao/module_ri/rpa_lri_sym.cpp b/source/source_lcao/module_ri/rpa_lri_sym.cpp
deleted file mode 100644
index fe6694b56ca..00000000000
--- a/source/source_lcao/module_ri/rpa_lri_sym.cpp
+++ /dev/null
@@ -1,69 +0,0 @@
-#include "rpa_lri_detail.h"
-#include "source_cell/klist.h"
-#include "source_cell/unitcell.h"
-
-#include <cmath>
-#include <iomanip>
-#include <ostream>
-#include <stdexcept>
-#include <vector>
-
-namespace
-{
-bool preserves_mesh(const ModuleBase::Matrix3& kg, const K_Vectors& kv, const UnitCell& cell)
-{
-    const double rows[3][3] = {{kg.e11, kg.e12, kg.e13}, {kg.e21, kg.e22, kg.e23}, {kg.e31, kg.e32, kg.e33}};
-    // ABACUS and LibRPA rotate row vectors: k' = k * kg.  Therefore the
-    // coefficient taking source component j on an nmp[j] mesh to target
-    // component i on an nmp[i] mesh is kg[j][i] * nmp[i] / nmp[j].
-    for (int j = 0; j < 3; ++j)
-    {
-        if (kv.nmp[j] <= 0) throw std::runtime_error("Spatial symmetry requires a uniform k-grid");
-        for (int i = 0; i < 3; ++i)
-        {
-            const double step = rows[j][i] * kv.nmp[i] / kv.nmp[j];
-            if (std::abs(step - std::round(step)) > 1e-8) return false;
-        }
-    }
-    // Check the actual mesh offset as well as its spacing.
-    const auto origin = kv.kvec_c_full.at(0) * cell.G.Inverse();
-    const auto delta = origin * kg - origin;
-    const double shifts[3] = {delta.x * kv.nmp[0], delta.y * kv.nmp[1], delta.z * kv.nmp[2]};
-    for (const double shift : shifts)
-    {
-        if (std::abs(shift - std::round(shift)) > 1e-8) return false;
-    }
-    return true;
-}
-}
-
-void RpaLriDetail::write_stru_sym(std::ostream& output, const UnitCell& cell, const K_Vectors& kv)
-{
-    const auto& symm = cell.symm;
-    std::vector<int> operations;
-    for (int isym = 0; isym < symm.nrotk; ++isym)
-    {
-        // Use the same reciprocal-space matrices that ABACUS used to reduce
-        // this mesh.  Reconstructing them from gmatrix can change the basis
-        // convention for non-orthogonal primitive cells.
-        if (preserves_mesh(symm.kgmatrix[isym], kv, cell)) operations.push_back(isym);
-    }
-    if (operations.empty()) throw std::runtime_error("Missing mesh-preserving spatial symmetry operations");
-    // Use the subgroup that preserves this mesh. The full crystal group may
-    // rotate an anisotropic mesh outside itself and cannot accelerate its RPA.
-    output << operations.size() << " row\n";
-    for (const int isym : operations)
-    {
-        const auto& r = symm.gmatrix[isym];
-        const double entries[9] = {r.e11, r.e12, r.e13, r.e21, r.e22, r.e23, r.e31, r.e32, r.e33};
-        for (const double entry : entries)
-        {
-            if (!std::isfinite(entry) || std::abs(entry - std::round(entry)) > 1e-8)
-                throw std::runtime_error("Nonintegral fractional symmetry rotation");
-            output << std::setw(4) << static_cast<int>(std::round(entry));
-        }
-        const auto& t = symm.gtrans[isym];
-        output << std::scientific << std::setprecision(15) << std::setw(24) << t.x
-               << std::setw(24) << t.y << std::setw(24) << t.z << '\n';
-    }
-}
diff --git a/source/source_lcao/module_ri/unittests/CMakeLists.txt b/source/source_lcao/module_ri/unittests/CMakeLists.txt
index 5635d5d68a1..d20a6ddb915 100644
--- a/source/source_lcao/module_ri/unittests/CMakeLists.txt
+++ b/source/source_lcao/module_ri/unittests/CMakeLists.txt
@@ -16,15 +16,5 @@ AddTest(
   LIBS parameter
   SOURCES test_abfs_vector3_order.cpp
 )
-AddTest(
-  TARGET MODULE_RI_rpa_lri_detail_test
-  LIBS parameter base device
-  SOURCES test_rpa_lri_detail.cpp ../rpa_lri_detail.cpp ../rpa_lri_coeff_v1.cpp ../rpa_lri_bz_io.cpp
-)
-AddTest(
-  TARGET MODULE_RI_rpa_lri_band_io_test
-  LIBS parameter base device
-  SOURCES test_rpa_lri_band_io.cpp ../rpa_lri_band_io.cpp
-)
 file(COPY "${CMAKE_CURRENT_SOURCE_DIR}/support"
      DESTINATION "${CMAKE_CURRENT_BINARY_DIR}")
diff --git a/source/source_lcao/module_rpa/CMakeLists.txt b/source/source_lcao/module_rpa/CMakeLists.txt
new file mode 100644
index 00000000000..9969d4431dd
--- /dev/null
+++ b/source/source_lcao/module_rpa/CMakeLists.txt
@@ -0,0 +1,42 @@
+if(NOT ENABLE_LIBRI OR NOT ENABLE_LCAO)
+  return()
+endif()
+
+add_subdirectory(module_librpa_io)
+
+add_library(rpa OBJECT
+  rpa_lri.cpp
+  rpa_lri_detail.cpp
+  rpa_lri_band_io.cpp
+  rpa_lri_coeff_v1.cpp
+  rpa_lri_bz_io.cpp
+  rpa_lri_mpi_layout.cpp
+  rpa_lri_wfc_mpi.cpp
+  rpa_lri_coulomb.cpp
+  rpa_lri_coulomb_io.cpp
+  rpa_lri_overlap.cpp
+  rpa_lri_overlap_io.cpp
+  rpa_lri_output.cpp
+  rpa_lri_bz.cpp
+  rpa_lri_sym.cpp
+  rpa_lri_coeff_io.cpp
+  rpa_lri_velocity.cpp
+  rpa_lri_basis.cpp
+  rpa_lri_eigenvector_v1.cpp
+  rpa_lri_eigenvector.cpp
+)
+target_include_directories(rpa PRIVATE
+  ${ABACUS_SOURCE_DIR}
+  ${ABACUS_SOURCE_DIR}/source_lcao/module_ri
+  ${CMAKE_CURRENT_SOURCE_DIR}
+  ${CMAKE_CURRENT_SOURCE_DIR}/module_librpa_io
+)
+target_link_libraries(rpa PRIVATE ri librpa_io)
+
+if(BUILD_TESTING AND ENABLE_MPI)
+  add_subdirectory(unittests)
+endif()
+
+if(ENABLE_COVERAGE)
+  add_coverage(rpa)
+endif()
diff --git a/source/source_lcao/module_rpa/module_librpa_io/CMakeLists.txt b/source/source_lcao/module_rpa/module_librpa_io/CMakeLists.txt
new file mode 100644
index 00000000000..202a7be663b
--- /dev/null
+++ b/source/source_lcao/module_rpa/module_librpa_io/CMakeLists.txt
@@ -0,0 +1,13 @@
+add_library(librpa_io OBJECT
+  librpa_bz_sampling.cpp
+  librpa_stru_symmetry.cpp
+  librpa_stru_units.cpp
+)
+target_include_directories(librpa_io PUBLIC
+  ${ABACUS_SOURCE_DIR}
+  ${CMAKE_CURRENT_SOURCE_DIR}
+)
+
+if(BUILD_TESTING AND ENABLE_MPI)
+  add_subdirectory(unittests)
+endif()
diff --git a/source/source_lcao/module_rpa/module_librpa_io/librpa_bz_sampling.cpp b/source/source_lcao/module_rpa/module_librpa_io/librpa_bz_sampling.cpp
new file mode 100644
index 00000000000..9e8eea6da14
--- /dev/null
+++ b/source/source_lcao/module_rpa/module_librpa_io/librpa_bz_sampling.cpp
@@ -0,0 +1,18 @@
+#include "librpa_bz_sampling.h"
+
+#include "source_base/constants.h"
+
+#include <cmath>
+#include <stdexcept>
+
+namespace RpaLriDetail
+{
+ModuleBase::Vector3<double> librpa_cartesian_kvec(const ModuleBase::Vector3<double>& kvec_c, const double lat0_bohr)
+{
+    if (!std::isfinite(lat0_bohr) || lat0_bohr <= 0.0)
+    {
+        throw std::invalid_argument("LibRPA BZ output requires a positive lattice constant.");
+    }
+    return kvec_c * (ModuleBase::TWO_PI / lat0_bohr);
+}
+} // namespace RpaLriDetail
diff --git a/source/source_lcao/module_rpa/module_librpa_io/librpa_bz_sampling.h b/source/source_lcao/module_rpa/module_librpa_io/librpa_bz_sampling.h
new file mode 100644
index 00000000000..6e8f4d6f4c9
--- /dev/null
+++ b/source/source_lcao/module_rpa/module_librpa_io/librpa_bz_sampling.h
@@ -0,0 +1,12 @@
+#ifndef LIBRPA_BZ_SAMPLING_H
+#define LIBRPA_BZ_SAMPLING_H
+
+#include "source_base/vector3.h"
+
+namespace RpaLriDetail
+{
+ModuleBase::Vector3<double> librpa_cartesian_kvec(const ModuleBase::Vector3<double>& kvec_c,
+                                                  double lat0_bohr);
+} // namespace RpaLriDetail
+
+#endif
diff --git a/source/source_lcao/module_rpa/module_librpa_io/librpa_stru_symmetry.cpp b/source/source_lcao/module_rpa/module_librpa_io/librpa_stru_symmetry.cpp
new file mode 100644
index 00000000000..db4d29a73e2
--- /dev/null
+++ b/source/source_lcao/module_rpa/module_librpa_io/librpa_stru_symmetry.cpp
@@ -0,0 +1,121 @@
+#include "librpa_stru_symmetry.h"
+
+#include "source_base/matrix3.h"
+#include "source_base/vector3.h"
+#include "source_cell/module_symmetry/symm_rot_spin.h"
+
+#include <cmath>
+#include <iomanip>
+#include <stdexcept>
+
+namespace RpaLriDetail
+{
+int checked_near_int(const double value, const std::string& context)
+{
+    const double rounded = std::round(value);
+    if (!std::isfinite(value) || std::abs(value - rounded) > 1e-8)
+    {
+        throw std::runtime_error(context + " is not close to an integer.");
+    }
+    return static_cast<int>(rounded);
+}
+
+LibRpaSymmetryOperation make_librpa_symmetry_operation(const ModuleBase::Matrix3& rotation,
+                                                       const ModuleBase::Vector3<double>& translation)
+{
+    LibRpaSymmetryOperation operation = {{rotation.e11,
+                                          rotation.e12,
+                                          rotation.e13,
+                                          rotation.e21,
+                                          rotation.e22,
+                                          rotation.e23,
+                                          rotation.e31,
+                                          rotation.e32,
+                                          rotation.e33},
+                                         {translation.x, translation.y, translation.z}};
+    return operation;
+}
+
+std::vector<LibRpaSpinSymmetryOperation> make_librpa_spin_symmetry_operations(
+    const ModuleBase::Matrix3& latvec,
+    const std::vector<ModuleBase::Matrix3>& unitary_rotations,
+    const std::vector<ModuleBase::Matrix3>& antiunitary_rotations)
+{
+    const auto spin_operation = [&latvec](const ModuleBase::Matrix3& rotation, const int antiunitary) {
+        // Convert the fractional row-vector rotation to Cartesian coordinates
+        // before deriving the SU(2) representation used by the spinor reader.
+        const ModuleBase::Matrix3 cartesian = latvec.Inverse() * rotation * latvec;
+        const auto spin_u = ModuleSymmetry::SpinRotation::so3_to_su2(cartesian);
+        LibRpaSpinSymmetryOperation operation;
+        operation.antiunitary = antiunitary;
+        for (int i = 0; i != 4; ++i)
+        {
+            operation.spin_u[static_cast<std::size_t>(2 * i)] = spin_u[i].real();
+            operation.spin_u[static_cast<std::size_t>(2 * i + 1)] = spin_u[i].imag();
+        }
+        return operation;
+    };
+
+    std::vector<LibRpaSpinSymmetryOperation> operations;
+    operations.reserve(unitary_rotations.size() + antiunitary_rotations.size());
+    for (const auto& rotation : unitary_rotations)
+    {
+        operations.push_back(spin_operation(rotation, 0));
+    }
+    for (const auto& rotation : antiunitary_rotations)
+    {
+        operations.push_back(spin_operation(rotation, 1));
+    }
+    return operations;
+}
+
+void write_librpa_symmetry_rows(std::ostream& output,
+                                const std::vector<LibRpaSymmetryOperation>& unitary,
+                                const std::vector<LibRpaSymmetryOperation>& antiunitary)
+{
+    const auto write_operation = [&output](const LibRpaSymmetryOperation& operation) {
+        output << std::setw(4) << checked_near_int(operation.rotation[0], "symmetry rotation e11")
+               << std::setw(4) << checked_near_int(operation.rotation[1], "symmetry rotation e12")
+               << std::setw(4) << checked_near_int(operation.rotation[2], "symmetry rotation e13")
+               << std::setw(4) << checked_near_int(operation.rotation[3], "symmetry rotation e21")
+               << std::setw(4) << checked_near_int(operation.rotation[4], "symmetry rotation e22")
+               << std::setw(4) << checked_near_int(operation.rotation[5], "symmetry rotation e23")
+               << std::setw(4) << checked_near_int(operation.rotation[6], "symmetry rotation e31")
+               << std::setw(4) << checked_near_int(operation.rotation[7], "symmetry rotation e32")
+               << std::setw(4) << checked_near_int(operation.rotation[8], "symmetry rotation e33")
+               << std::setw(24) << std::scientific << std::setprecision(15) << operation.translation[0]
+               << std::setw(24) << std::scientific << std::setprecision(15) << operation.translation[1]
+               << std::setw(24) << std::scientific << std::setprecision(15) << operation.translation[2] << '\n';
+    };
+
+    output << unitary.size() + antiunitary.size() << " row\n";
+    for (const auto& operation : unitary)
+    {
+        write_operation(operation);
+    }
+    for (const auto& operation : antiunitary)
+    {
+        write_operation(operation);
+    }
+}
+
+void write_librpa_spin_symmetry(std::ostream& output,
+                                const int grey_group,
+                                const int spin_source,
+                                const std::vector<LibRpaSpinSymmetryOperation>& operations)
+{
+    output << "spin_symmetry " << grey_group << ' ' << spin_source << '\n';
+    for (const auto& operation : operations)
+    {
+        output << operation.antiunitary;
+        if (spin_source == 1)
+        {
+            for (const double value : operation.spin_u)
+            {
+                output << std::setw(24) << std::scientific << std::setprecision(15) << value;
+            }
+        }
+        output << '\n';
+    }
+}
+} // namespace RpaLriDetail
diff --git a/source/source_lcao/module_rpa/module_librpa_io/librpa_stru_symmetry.h b/source/source_lcao/module_rpa/module_librpa_io/librpa_stru_symmetry.h
new file mode 100644
index 00000000000..f48f355ef48
--- /dev/null
+++ b/source/source_lcao/module_rpa/module_librpa_io/librpa_stru_symmetry.h
@@ -0,0 +1,50 @@
+#ifndef LIBRPA_STRU_SYMMETRY_H
+#define LIBRPA_STRU_SYMMETRY_H
+
+#include <array>
+#include <iosfwd>
+#include <string>
+#include <vector>
+
+namespace ModuleBase
+{
+class Matrix3;
+template <typename T>
+class Vector3;
+}
+
+namespace RpaLriDetail
+{
+struct LibRpaSymmetryOperation
+{
+    double rotation[9];
+    double translation[3];
+};
+
+struct LibRpaSpinSymmetryOperation
+{
+    int antiunitary = 0;
+    std::array<double, 8> spin_u{};
+};
+
+int checked_near_int(double value, const std::string& context);
+
+LibRpaSymmetryOperation make_librpa_symmetry_operation(const ModuleBase::Matrix3& rotation,
+                                                       const ModuleBase::Vector3<double>& translation);
+
+std::vector<LibRpaSpinSymmetryOperation> make_librpa_spin_symmetry_operations(
+    const ModuleBase::Matrix3& latvec,
+    const std::vector<ModuleBase::Matrix3>& unitary_rotations,
+    const std::vector<ModuleBase::Matrix3>& antiunitary_rotations);
+
+void write_librpa_symmetry_rows(std::ostream& output,
+                                const std::vector<LibRpaSymmetryOperation>& unitary,
+                                const std::vector<LibRpaSymmetryOperation>& antiunitary);
+
+void write_librpa_spin_symmetry(std::ostream& output,
+                                int grey_group,
+                                int spin_source,
+                                const std::vector<LibRpaSpinSymmetryOperation>& operations);
+} // namespace RpaLriDetail
+
+#endif
diff --git a/source/source_lcao/module_rpa/module_librpa_io/librpa_stru_units.cpp b/source/source_lcao/module_rpa/module_librpa_io/librpa_stru_units.cpp
new file mode 100644
index 00000000000..fc8cf6c2652
--- /dev/null
+++ b/source/source_lcao/module_rpa/module_librpa_io/librpa_stru_units.cpp
@@ -0,0 +1,18 @@
+#include "librpa_stru_units.h"
+
+#include "source_base/constants.h"
+
+#include <cmath>
+#include <stdexcept>
+
+namespace RpaLriDetail
+{
+LibRpaStruUnitScales librpa_stru_unit_scales(const double lat0_bohr)
+{
+    if (!std::isfinite(lat0_bohr) || lat0_bohr <= 0.0)
+    {
+        throw std::invalid_argument("LibRPA stru_out requires a positive finite lattice constant.");
+    }
+    return {lat0_bohr, ModuleBase::TWO_PI / lat0_bohr};
+}
+} // namespace RpaLriDetail
diff --git a/source/source_lcao/module_rpa/module_librpa_io/librpa_stru_units.h b/source/source_lcao/module_rpa/module_librpa_io/librpa_stru_units.h
new file mode 100644
index 00000000000..4d14e4d19cf
--- /dev/null
+++ b/source/source_lcao/module_rpa/module_librpa_io/librpa_stru_units.h
@@ -0,0 +1,15 @@
+#ifndef LIBRPA_STRU_UNITS_H
+#define LIBRPA_STRU_UNITS_H
+
+namespace RpaLriDetail
+{
+struct LibRpaStruUnitScales
+{
+    double real_space_bohr;
+    double reciprocal_space_bohr_inv;
+};
+
+LibRpaStruUnitScales librpa_stru_unit_scales(double lat0_bohr);
+} // namespace RpaLriDetail
+
+#endif
diff --git a/source/source_lcao/module_rpa/module_librpa_io/unittests/CMakeLists.txt b/source/source_lcao/module_rpa/module_librpa_io/unittests/CMakeLists.txt
new file mode 100644
index 00000000000..bb351373876
--- /dev/null
+++ b/source/source_lcao/module_rpa/module_librpa_io/unittests/CMakeLists.txt
@@ -0,0 +1,8 @@
+AddTest(
+  TARGET MODULE_RPA_librpa_stru_symmetry_test
+  LIBS base device symmetry
+  SOURCES test_librpa_stru_symmetry.cpp
+          ../librpa_bz_sampling.cpp
+          ../librpa_stru_symmetry.cpp
+          ../librpa_stru_units.cpp
+)
diff --git a/source/source_lcao/module_rpa/module_librpa_io/unittests/test_librpa_stru_symmetry.cpp b/source/source_lcao/module_rpa/module_librpa_io/unittests/test_librpa_stru_symmetry.cpp
new file mode 100644
index 00000000000..351845b22ed
--- /dev/null
+++ b/source/source_lcao/module_rpa/module_librpa_io/unittests/test_librpa_stru_symmetry.cpp
@@ -0,0 +1,132 @@
+#include "../librpa_stru_symmetry.h"
+#include "../librpa_stru_units.h"
+#include "../librpa_bz_sampling.h"
+
+#include "source_base/constants.h"
+#include "source_base/matrix3.h"
+#include "source_base/vector3.h"
+#include "source_cell/module_symmetry/symm_rot_spin.h"
+
+#include <cmath>
+#include <gtest/gtest.h>
+#include <limits>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+namespace
+{
+TEST(RpaLriSymmetry, WritesUnitaryRowsBeforeAntiunitaryRows)
+{
+    const std::vector<RpaLriDetail::LibRpaSymmetryOperation> unitary = {
+        RpaLriDetail::make_librpa_symmetry_operation(ModuleBase::Matrix3(), ModuleBase::Vector3<double>())};
+    const std::vector<RpaLriDetail::LibRpaSymmetryOperation> antiunitary = {
+        RpaLriDetail::make_librpa_symmetry_operation(
+            ModuleBase::Matrix3(-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0),
+            ModuleBase::Vector3<double>(0.5, 0.0, 0.0))};
+
+    std::ostringstream output;
+    RpaLriDetail::write_librpa_symmetry_rows(output, unitary, antiunitary);
+    const std::string text = output.str();
+    std::istringstream lines(text);
+    std::string line;
+    std::vector<std::string> rows;
+    while (std::getline(lines, line))
+    {
+        rows.push_back(line);
+    }
+
+    ASSERT_EQ(rows.size(), 3u);
+    EXPECT_EQ(rows[0], "2 row");
+    std::istringstream unitary_row(rows[1]);
+    std::istringstream antiunitary_row(rows[2]);
+    int unitary_rotation[9] = {};
+    int antiunitary_rotation[9] = {};
+    for (int i = 0; i != 9; ++i)
+    {
+        unitary_row >> unitary_rotation[i];
+        antiunitary_row >> antiunitary_rotation[i];
+    }
+    EXPECT_EQ(antiunitary_rotation[0], -1);
+    for (int i = 0; i != 9; ++i)
+    {
+        const int expected = i == 0 || i == 4 || i == 8 ? 1 : 0;
+        EXPECT_EQ(unitary_rotation[i], expected);
+    }
+    EXPECT_EQ(antiunitary_rotation[4], -1);
+    EXPECT_EQ(antiunitary_rotation[8], 1);
+    EXPECT_NE(rows[2].find("5.000000000000000e-01"), std::string::npos);
+}
+
+TEST(RpaLriSymmetry, WritesSpinRowsWithAntiunitaryFlag)
+{
+    const std::vector<ModuleBase::Matrix3> unitary = {ModuleBase::Matrix3()};
+    const std::vector<ModuleBase::Matrix3> antiunitary = {
+        ModuleBase::Matrix3(0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0)};
+    const auto operations
+        = RpaLriDetail::make_librpa_spin_symmetry_operations(ModuleBase::Matrix3(), unitary, antiunitary);
+
+    ASSERT_EQ(operations.size(), 2u);
+    EXPECT_EQ(operations[0].antiunitary, 0);
+    EXPECT_EQ(operations[1].antiunitary, 1);
+    EXPECT_NEAR(operations[0].spin_u[0], 1.0, 1e-12);
+    EXPECT_NEAR(operations[0].spin_u[7], 0.0, 1e-12);
+
+    std::ostringstream output;
+    RpaLriDetail::write_librpa_spin_symmetry(output, 0, 1, operations);
+    const std::string text = output.str();
+    EXPECT_NE(text.find("spin_symmetry 0 1\n"), std::string::npos);
+    EXPECT_NE(text.find("\n1 "), std::string::npos);
+}
+
+TEST(RpaLriSymmetry, UsesPhysicalLatticeConstantForBothUnits)
+{
+    const auto scales = RpaLriDetail::librpa_stru_unit_scales(5.3);
+    EXPECT_DOUBLE_EQ(scales.real_space_bohr, 5.3);
+    EXPECT_DOUBLE_EQ(scales.reciprocal_space_bohr_inv, ModuleBase::TWO_PI / 5.3);
+    EXPECT_NEAR(scales.real_space_bohr * scales.reciprocal_space_bohr_inv,
+                ModuleBase::TWO_PI,
+                4.0 * std::numeric_limits<double>::epsilon());
+    EXPECT_THROW(RpaLriDetail::librpa_stru_unit_scales(0.0), std::invalid_argument);
+}
+
+TEST(RpaLriSymmetry, ConvertsCartesianKPointWithPhysicalScale)
+{
+    const auto cartesian = RpaLriDetail::librpa_cartesian_kvec(ModuleBase::Vector3<double>(1.0, -2.0, 0.5), 2.0);
+    EXPECT_DOUBLE_EQ(cartesian.x, ModuleBase::PI);
+    EXPECT_DOUBLE_EQ(cartesian.y, -2.0 * ModuleBase::PI);
+    EXPECT_DOUBLE_EQ(cartesian.z, 0.5 * ModuleBase::PI);
+    EXPECT_THROW(RpaLriDetail::librpa_cartesian_kvec(ModuleBase::Vector3<double>(), -1.0), std::invalid_argument);
+}
+
+TEST(RpaLriSymmetry, ConvertsNonorthogonalLatticeRotationBeforeSpinMapping)
+{
+    const double half = 0.5;
+    const double height = std::sqrt(3.0) * half;
+    const ModuleBase::Matrix3 latvec(1.0, 0.0, 0.0, half, height, 0.0, 0.0, 0.0, 1.0);
+    // This fractional operation is a 60-degree rotation of the triangular lattice.
+    const ModuleBase::Matrix3 fractional_rotation(0.0, 1.0, 0.0, -1.0, 1.0, 0.0, 0.0, 0.0, 1.0);
+    const auto operations = RpaLriDetail::make_librpa_spin_symmetry_operations(
+        latvec, std::vector<ModuleBase::Matrix3>{fractional_rotation}, {});
+
+    ASSERT_EQ(operations.size(), 1u);
+    const double half_angle_cos = std::sqrt(3.0) * half;
+    EXPECT_NEAR(operations[0].spin_u[0], half_angle_cos, 1e-12);
+    EXPECT_NEAR(operations[0].spin_u[1], -half, 1e-12);
+    EXPECT_NEAR(operations[0].spin_u[2], 0.0, 1e-12);
+    EXPECT_NEAR(operations[0].spin_u[3], 0.0, 1e-12);
+    EXPECT_NEAR(operations[0].spin_u[4], 0.0, 1e-12);
+    EXPECT_NEAR(operations[0].spin_u[5], 0.0, 1e-12);
+    EXPECT_NEAR(operations[0].spin_u[6], half_angle_cos, 1e-12);
+    EXPECT_NEAR(operations[0].spin_u[7], half, 1e-12);
+
+    const ModuleBase::Matrix3 cartesian = latvec.Inverse() * fractional_rotation * latvec;
+    const auto expected = ModuleSymmetry::SpinRotation::so3_to_su2(cartesian);
+    for (int i = 0; i != 4; ++i)
+    {
+        EXPECT_NEAR(operations[0].spin_u[static_cast<std::size_t>(2 * i)], expected[i].real(), 1e-12);
+        EXPECT_NEAR(operations[0].spin_u[static_cast<std::size_t>(2 * i + 1)], expected[i].imag(), 1e-12);
+    }
+}
+} // namespace
diff --git a/source/source_lcao/module_ri/rpa_lri.cpp b/source/source_lcao/module_rpa/rpa_lri.cpp
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri.cpp
rename to source/source_lcao/module_rpa/rpa_lri.cpp
diff --git a/source/source_lcao/module_ri/rpa_lri.h b/source/source_lcao/module_rpa/rpa_lri.h
similarity index 99%
rename from source/source_lcao/module_ri/rpa_lri.h
rename to source/source_lcao/module_rpa/rpa_lri.h
index 5ac44987ca3..e9a98c0bf53 100644
--- a/source/source_lcao/module_ri/rpa_lri.h
+++ b/source/source_lcao/module_rpa/rpa_lri.h
@@ -6,7 +6,7 @@
 #ifndef RPA_LRI_H
 #define RPA_LRI_H
 
-#include "lri_cv.h"
+#include "source_lcao/module_ri/lri_cv.h"
 #include "source_hamilt/module_xc/exx_info_ri.h"
 #include "source_psi/psi.h"
 // #include "module_xc/exx_info.h"
diff --git a/source/source_lcao/module_ri/rpa_lri_band_io.cpp b/source/source_lcao/module_rpa/rpa_lri_band_io.cpp
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_band_io.cpp
rename to source/source_lcao/module_rpa/rpa_lri_band_io.cpp
diff --git a/source/source_lcao/module_ri/rpa_lri_band_io.h b/source/source_lcao/module_rpa/rpa_lri_band_io.h
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_band_io.h
rename to source/source_lcao/module_rpa/rpa_lri_band_io.h
diff --git a/source/source_lcao/module_ri/rpa_lri_basis.cpp b/source/source_lcao/module_rpa/rpa_lri_basis.cpp
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_basis.cpp
rename to source/source_lcao/module_rpa/rpa_lri_basis.cpp
diff --git a/source/source_lcao/module_ri/rpa_lri_bz.cpp b/source/source_lcao/module_rpa/rpa_lri_bz.cpp
similarity index 87%
rename from source/source_lcao/module_ri/rpa_lri_bz.cpp
rename to source/source_lcao/module_rpa/rpa_lri_bz.cpp
index 9b255e19557..f76aacc5620 100644
--- a/source/source_lcao/module_ri/rpa_lri_bz.cpp
+++ b/source/source_lcao/module_rpa/rpa_lri_bz.cpp
@@ -1,6 +1,7 @@
 #include "rpa_lri.h"
 #include "rpa_lri_bz_io.h"
 #include "rpa_lri_detail.h"
+#include "librpa_bz_sampling.h"
 #include "source_base/global_function.h"
 #include "source_cell/klist.h"
 #include "source_cell/unitcell.h"
@@ -26,7 +27,7 @@ try
     data.reduced.resize(nks_tot);
     for (int ik = 0; ik < nks_tot; ++ik)
     {
-        const ModuleBase::Vector3<double> kcart = p_kv->kvec_c[ik] * (ModuleBase::TWO_PI / ucell.lat0);
+        const ModuleBase::Vector3<double> kcart = RpaLriDetail::librpa_cartesian_kvec(p_kv->kvec_c[ik], ucell.lat0);
         data.reduced[ik].weight = p_kv->wk[ik];
         data.reduced[ik].direct = {{p_kv->kvec_d[ik].x, p_kv->kvec_d[ik].y, p_kv->kvec_d[ik].z}};
         data.reduced[ik].cartesian = {{kcart.x, kcart.y, kcart.z}};
@@ -34,7 +35,7 @@ try
     data.full_cartesian.resize(p_kv->kvec_c_full.size());
     for (std::size_t ik = 0; ik < p_kv->kvec_c_full.size(); ++ik)
     {
-        const auto kcart = p_kv->kvec_c_full[ik] * (ModuleBase::TWO_PI / ucell.lat0);
+        const auto kcart = RpaLriDetail::librpa_cartesian_kvec(p_kv->kvec_c_full[ik], ucell.lat0);
         data.full_cartesian[ik] = {{kcart.x, kcart.y, kcart.z}};
     }
     data.full_map = p_kv->ibz_index;
diff --git a/source/source_lcao/module_ri/rpa_lri_bz_io.cpp b/source/source_lcao/module_rpa/rpa_lri_bz_io.cpp
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_bz_io.cpp
rename to source/source_lcao/module_rpa/rpa_lri_bz_io.cpp
diff --git a/source/source_lcao/module_ri/rpa_lri_bz_io.h b/source/source_lcao/module_rpa/rpa_lri_bz_io.h
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_bz_io.h
rename to source/source_lcao/module_rpa/rpa_lri_bz_io.h
diff --git a/source/source_lcao/module_ri/rpa_lri_coeff_io.cpp b/source/source_lcao/module_rpa/rpa_lri_coeff_io.cpp
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_coeff_io.cpp
rename to source/source_lcao/module_rpa/rpa_lri_coeff_io.cpp
diff --git a/source/source_lcao/module_ri/rpa_lri_coeff_v1.cpp b/source/source_lcao/module_rpa/rpa_lri_coeff_v1.cpp
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_coeff_v1.cpp
rename to source/source_lcao/module_rpa/rpa_lri_coeff_v1.cpp
diff --git a/source/source_lcao/module_ri/rpa_lri_coeff_v1.h b/source/source_lcao/module_rpa/rpa_lri_coeff_v1.h
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_coeff_v1.h
rename to source/source_lcao/module_rpa/rpa_lri_coeff_v1.h
diff --git a/source/source_lcao/module_ri/rpa_lri_coulomb.cpp b/source/source_lcao/module_rpa/rpa_lri_coulomb.cpp
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_coulomb.cpp
rename to source/source_lcao/module_rpa/rpa_lri_coulomb.cpp
diff --git a/source/source_lcao/module_ri/rpa_lri_coulomb_io.cpp b/source/source_lcao/module_rpa/rpa_lri_coulomb_io.cpp
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_coulomb_io.cpp
rename to source/source_lcao/module_rpa/rpa_lri_coulomb_io.cpp
diff --git a/source/source_lcao/module_ri/rpa_lri_detail.cpp b/source/source_lcao/module_rpa/rpa_lri_detail.cpp
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_detail.cpp
rename to source/source_lcao/module_rpa/rpa_lri_detail.cpp
diff --git a/source/source_lcao/module_ri/rpa_lri_detail.h b/source/source_lcao/module_rpa/rpa_lri_detail.h
similarity index 99%
rename from source/source_lcao/module_ri/rpa_lri_detail.h
rename to source/source_lcao/module_rpa/rpa_lri_detail.h
index d2b157c0120..6d40779ab44 100644
--- a/source/source_lcao/module_ri/rpa_lri_detail.h
+++ b/source/source_lcao/module_rpa/rpa_lri_detail.h
@@ -34,7 +34,7 @@ static_assert(sizeof(std::complex<double>) == 2 * sizeof(double),
               "LibRPA v1 binary output expects complex<double> as two doubles.");
 
 void trim_malloc_cache();
-void write_stru_sym(std::ostream& output, const UnitCell& cell, const K_Vectors& kv);
+void write_stru_sym(std::ostream& output, const UnitCell& cell, const K_Vectors& kv, const int nspin);
 
 bool debug_dump_exx_ao_enabled();
 
diff --git a/source/source_lcao/module_ri/rpa_lri_eigenvector.cpp b/source/source_lcao/module_rpa/rpa_lri_eigenvector.cpp
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_eigenvector.cpp
rename to source/source_lcao/module_rpa/rpa_lri_eigenvector.cpp
diff --git a/source/source_lcao/module_ri/rpa_lri_eigenvector_v1.cpp b/source/source_lcao/module_rpa/rpa_lri_eigenvector_v1.cpp
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_eigenvector_v1.cpp
rename to source/source_lcao/module_rpa/rpa_lri_eigenvector_v1.cpp
diff --git a/source/source_lcao/module_ri/rpa_lri_mpi_layout.cpp b/source/source_lcao/module_rpa/rpa_lri_mpi_layout.cpp
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_mpi_layout.cpp
rename to source/source_lcao/module_rpa/rpa_lri_mpi_layout.cpp
diff --git a/source/source_lcao/module_ri/rpa_lri_output.cpp b/source/source_lcao/module_rpa/rpa_lri_output.cpp
similarity index 93%
rename from source/source_lcao/module_ri/rpa_lri_output.cpp
rename to source/source_lcao/module_rpa/rpa_lri_output.cpp
index 5b3289ba16d..4dc970cc3d1 100644
--- a/source/source_lcao/module_ri/rpa_lri_output.cpp
+++ b/source/source_lcao/module_rpa/rpa_lri_output.cpp
@@ -2,6 +2,7 @@
 #include "rpa_lri.h"
 #include "rpa_lri_band_io.h"
 #include "rpa_lri_detail.h"
+#include "librpa_stru_units.h"
 #include "source_base/global_function.h"
 #include "source_basis/module_ao/elem_basis_idx_orb.h"
 #include "source_estate/elecstate_lcao.h"
@@ -34,8 +35,9 @@ try
     }
     ModuleBase::TITLE("DFT_RPA_interface", "out_struc");
     const int nks_tot = this->runtime.input.nspin == 2 ? (int)p_kv->get_nks() / 2 : p_kv->get_nks();
-    const ModuleBase::Matrix3 lat = ucell.latvec * ucell.lat0;                     // in unit of Bohr
-    const ModuleBase::Matrix3 G_RPA = ucell.G * (ModuleBase::TWO_PI / ucell.lat0); // in unit of 1/Bohr
+    const auto unit_scales = RpaLriDetail::librpa_stru_unit_scales(ucell.lat0);
+    const ModuleBase::Matrix3 lat = ucell.latvec * unit_scales.real_space_bohr;
+    const ModuleBase::Matrix3 G_RPA = ucell.G * unit_scales.reciprocal_space_bohr_inv;
     std::ofstream ofs;
     ofs.open(this->runtime.input.rpa_outdir + "stru_out.txt", std::ios::out);
     if (!ofs.good()) throw std::runtime_error("Failed to open stru_out.txt");
@@ -53,7 +55,7 @@ try
     {
         for (int ia = 0; ia < ucell.atoms[it].na; ia++)
         {
-            const ModuleBase::Vector3<double> position = ucell.atoms[it].tau[ia] * ucell.lat0; // in unit of Bohr
+            const ModuleBase::Vector3<double> position = ucell.atoms[it].tau[ia] * unit_scales.real_space_bohr;
             ofs << std::setw(15) << position.x << std::setw(15) << position.y << std::setw(15) << position.z
                 << std::setw(15) << it + 1 << std::endl;
         }
@@ -84,10 +86,9 @@ try
     // the operations in the input cell's fractional row-vector convention.
     // This block belongs only to the reader-v1 structure contract; v0 keeps
     // the historical stru_out layout unchanged.
-    if (this->runtime.input.out_librpa_ver == 1 && this->runtime.input.symmetry == "1"
-        && this->runtime.input.nspin != 4)
+    if (this->runtime.input.out_librpa_ver == 1 && this->runtime.input.symmetry == "1")
     {
-        RpaLriDetail::write_stru_sym(ofs, ucell, *p_kv);
+        RpaLriDetail::write_stru_sym(ofs, ucell, *p_kv, this->runtime.input.nspin);
     }
     ofs.close();
     if (!ofs.good()) throw std::runtime_error("Failed to write stru_out.txt");
diff --git a/source/source_lcao/module_ri/rpa_lri_overlap.cpp b/source/source_lcao/module_rpa/rpa_lri_overlap.cpp
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_overlap.cpp
rename to source/source_lcao/module_rpa/rpa_lri_overlap.cpp
diff --git a/source/source_lcao/module_ri/rpa_lri_overlap_io.cpp b/source/source_lcao/module_rpa/rpa_lri_overlap_io.cpp
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_overlap_io.cpp
rename to source/source_lcao/module_rpa/rpa_lri_overlap_io.cpp
diff --git a/source/source_lcao/module_rpa/rpa_lri_sym.cpp b/source/source_lcao/module_rpa/rpa_lri_sym.cpp
new file mode 100644
index 00000000000..018cd428ec2
--- /dev/null
+++ b/source/source_lcao/module_rpa/rpa_lri_sym.cpp
@@ -0,0 +1,94 @@
+#include "rpa_lri_detail.h"
+#include "librpa_stru_symmetry.h"
+#include "source_cell/klist.h"
+#include "source_cell/unitcell.h"
+
+#include <cmath>
+#include <iomanip>
+#include <ostream>
+#include <stdexcept>
+#include <vector>
+
+namespace
+{
+bool preserves_mesh(const ModuleBase::Matrix3& kg, const K_Vectors& kv, const UnitCell& cell)
+{
+    const double rows[3][3] = {{kg.e11, kg.e12, kg.e13}, {kg.e21, kg.e22, kg.e23}, {kg.e31, kg.e32, kg.e33}};
+    // ABACUS and LibRPA rotate row vectors: k' = k * kg.  Therefore the
+    // coefficient taking source component j on an nmp[j] mesh to target
+    // component i on an nmp[i] mesh is kg[j][i] * nmp[i] / nmp[j].
+    for (int j = 0; j < 3; ++j)
+    {
+        if (kv.nmp[j] <= 0) throw std::runtime_error("Spatial symmetry requires a uniform k-grid");
+        for (int i = 0; i < 3; ++i)
+        {
+            const double step = rows[j][i] * kv.nmp[i] / kv.nmp[j];
+            if (std::abs(step - std::round(step)) > 1e-8) return false;
+        }
+    }
+    // Check the actual mesh offset as well as its spacing.
+    const auto origin = kv.kvec_c_full.at(0) * cell.G.Inverse();
+    const auto delta = origin * kg - origin;
+    const double shifts[3] = {delta.x * kv.nmp[0], delta.y * kv.nmp[1], delta.z * kv.nmp[2]};
+    for (const double shift : shifts)
+    {
+        if (std::abs(shift - std::round(shift)) > 1e-8) return false;
+    }
+    return true;
+}
+
+ModuleBase::Matrix3 anti_k_action(const ModuleBase::Matrix3& spatial)
+{
+    ModuleBase::Matrix3 action = spatial;
+    action *= -1.0;
+    return action;
+}
+}
+
+void RpaLriDetail::write_stru_sym(std::ostream& output, const UnitCell& cell, const K_Vectors& kv, const int nspin)
+{
+    const auto& symm = cell.symm;
+    std::vector<LibRpaSymmetryOperation> unitary_operations;
+    std::vector<LibRpaSymmetryOperation> antiunitary_operations;
+    std::vector<ModuleBase::Matrix3> unitary_rotations;
+    std::vector<ModuleBase::Matrix3> antiunitary_rotations;
+    for (int isym = 0; isym < symm.nrotk; ++isym)
+    {
+        // Use the same reciprocal-space matrices that ABACUS used to reduce
+        // this mesh.  Reconstructing them from gmatrix can change the basis
+        // convention for non-orthogonal primitive cells.
+        if (preserves_mesh(symm.kgmatrix[isym], kv, cell))
+        {
+            unitary_operations.push_back(make_librpa_symmetry_operation(symm.gmatrix[isym], symm.gtrans[isym]));
+            unitary_rotations.push_back(symm.gmatrix[isym]);
+        }
+    }
+    if (unitary_operations.empty())
+    {
+        throw std::runtime_error("Missing mesh-preserving spatial symmetry operations");
+    }
+
+    if (nspin == 4 && symm.magnetic_nspin4)
+    {
+        for (int isym = 0; isym < symm.nrotk_anti; ++isym)
+        {
+            // The exported row stores the spatial part g. The k-point action
+            // of Theta*g is -k*g, so the mesh test must include time reversal.
+            if (preserves_mesh(anti_k_action(symm.kgmatrix_anti[isym]), kv, cell))
+            {
+                antiunitary_operations.push_back(
+                    make_librpa_symmetry_operation(symm.gmatrix_anti[isym], symm.gtrans_anti[isym]));
+                antiunitary_rotations.push_back(symm.gmatrix_anti[isym]);
+            }
+        }
+    }
+
+    write_librpa_symmetry_rows(output, unitary_operations, antiunitary_operations);
+    if (nspin == 4)
+    {
+        const auto spin_operations
+            = make_librpa_spin_symmetry_operations(cell.latvec, unitary_rotations, antiunitary_rotations);
+        const int grey_group = symm.magnetic_nspin4 ? 0 : 1;
+        write_librpa_spin_symmetry(output, grey_group, 1, spin_operations);
+    }
+}
diff --git a/source/source_lcao/module_ri/rpa_lri_velocity.cpp b/source/source_lcao/module_rpa/rpa_lri_velocity.cpp
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_velocity.cpp
rename to source/source_lcao/module_rpa/rpa_lri_velocity.cpp
diff --git a/source/source_lcao/module_ri/rpa_lri_wfc_mpi.cpp b/source/source_lcao/module_rpa/rpa_lri_wfc_mpi.cpp
similarity index 100%
rename from source/source_lcao/module_ri/rpa_lri_wfc_mpi.cpp
rename to source/source_lcao/module_rpa/rpa_lri_wfc_mpi.cpp
diff --git a/source/source_lcao/module_rpa/unittests/CMakeLists.txt b/source/source_lcao/module_rpa/unittests/CMakeLists.txt
new file mode 100644
index 00000000000..cc94ec19e45
--- /dev/null
+++ b/source/source_lcao/module_rpa/unittests/CMakeLists.txt
@@ -0,0 +1,14 @@
+AddTest(
+  TARGET MODULE_RPA_rpa_lri_detail_test
+  LIBS parameter base device
+  SOURCES test_rpa_lri_detail.cpp
+          ../rpa_lri_detail.cpp
+          ../rpa_lri_coeff_v1.cpp
+          ../rpa_lri_bz_io.cpp
+          ../module_librpa_io/librpa_bz_sampling.cpp
+)
+AddTest(
+  TARGET MODULE_RPA_rpa_lri_band_io_test
+  LIBS parameter base device
+  SOURCES test_rpa_lri_band_io.cpp ../rpa_lri_band_io.cpp
+)
diff --git a/source/source_lcao/module_ri/unittests/test_rpa_lri_band_io.cpp b/source/source_lcao/module_rpa/unittests/test_rpa_lri_band_io.cpp
similarity index 100%
rename from source/source_lcao/module_ri/unittests/test_rpa_lri_band_io.cpp
rename to source/source_lcao/module_rpa/unittests/test_rpa_lri_band_io.cpp
diff --git a/source/source_lcao/module_ri/unittests/test_rpa_lri_detail.cpp b/source/source_lcao/module_rpa/unittests/test_rpa_lri_detail.cpp
similarity index 95%
rename from source/source_lcao/module_ri/unittests/test_rpa_lri_detail.cpp
rename to source/source_lcao/module_rpa/unittests/test_rpa_lri_detail.cpp
index bb2894c22c1..2e5a06efa4e 100644
--- a/source/source_lcao/module_ri/unittests/test_rpa_lri_detail.cpp
+++ b/source/source_lcao/module_rpa/unittests/test_rpa_lri_detail.cpp
@@ -1,7 +1,9 @@
 #include "../rpa_lri_detail.h"
+#include "source_base/constants.h"
 
 #include "../rpa_lri_bz_io.h"
 #include "../rpa_lri_coeff_v1.h"
+#include "../module_librpa_io/librpa_bz_sampling.h"
 
 #include <RI/global/Tensor.h>
 #include <array>
@@ -35,6 +37,17 @@ using RpaLriDetail::write_bz_sample;
 using RpaLriDetail::write_lri_coef_v1;
 }
 
+TEST(RpaLriDetail, BzSamplingConvertsCartesianVectorsToInverseBohr)
+{
+    const auto cartesian
+        = RpaLriDetail::librpa_cartesian_kvec(ModuleBase::Vector3<double>(0.25, -0.5, 0.0), 2.0);
+    const double scale = ModuleBase::TWO_PI / 2.0;
+    EXPECT_DOUBLE_EQ(cartesian.x, 0.25 * scale);
+    EXPECT_DOUBLE_EQ(cartesian.y, -0.5 * scale);
+    EXPECT_DOUBLE_EQ(cartesian.z, 0.0);
+    EXPECT_THROW(RpaLriDetail::librpa_cartesian_kvec(ModuleBase::Vector3<double>(), 0.0), std::invalid_argument);
+}
+
 TEST(RpaLriDetail, CoulombAtomPairIndexUsesUpperTriangleOrder)
 {
     const std::size_t natoms = 4;
diff --git a/tests/08_RI/63_librpa_soc_symmetry/INPUT b/tests/08_RI/63_librpa_soc_symmetry/INPUT
new file mode 100644
index 00000000000..343a8eb615f
--- /dev/null
+++ b/tests/08_RI/63_librpa_soc_symmetry/INPUT
@@ -0,0 +1,30 @@
+INPUT_PARAMETERS
+suffix            autotest
+pseudo_dir        ../../PP_ORB
+orbital_dir       ../../PP_ORB
+calculation       scf
+ntype             1
+nbands            8
+gamma_only        0
+ecutwfc           20
+lspinorb          1
+nspin             4
+smearing_method   gaussian
+smearing_sigma    1e-4
+basis_type        lcao
+ks_solver         genelpa
+dft_functional    pbe
+rpa               1
+rpa_out_vel       1
+out_mat_xc        1
+out_wfc_lcao      1
+out_librpa_ver    1
+symmetry          1
+exx_separate_loop 1
+exx_pca_threshold 1e-3
+shrink_abfs_pca_thr 1e-1
+shrink_lu_inv_thr 1e-3
+exx_singularity_correction massidda
+exx_ccp_rmesh_times 1
+rpa_ccp_rmesh_times 1
+exx_real_number   1
diff --git a/tests/08_RI/63_librpa_soc_symmetry/KPT b/tests/08_RI/63_librpa_soc_symmetry/KPT
new file mode 100644
index 00000000000..a427751e525
--- /dev/null
+++ b/tests/08_RI/63_librpa_soc_symmetry/KPT
@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+1 1 2 0 0 0
diff --git a/tests/08_RI/63_librpa_soc_symmetry/README b/tests/08_RI/63_librpa_soc_symmetry/README
new file mode 100644
index 00000000000..ab723840b08
--- /dev/null
+++ b/tests/08_RI/63_librpa_soc_symmetry/README
@@ -0,0 +1,5 @@
+Check the reader-v1 structure export for a compact SOC calculation with
+symmetry enabled. The manifest verifies the spatial and spin symmetry tables;
+the test does not compare energies between symmetry settings.
+The symmetry runner also creates a shifted magnetic variant to exercise the
+antiunitary k-point mesh filter.
diff --git a/tests/08_RI/63_librpa_soc_symmetry/STRU b/tests/08_RI/63_librpa_soc_symmetry/STRU
new file mode 100644
index 00000000000..3d8eb4fdd15
--- /dev/null
+++ b/tests/08_RI/63_librpa_soc_symmetry/STRU
@@ -0,0 +1,22 @@
+ATOMIC_SPECIES
+H 1.008 H_ONCV_PBE_FR-1.0.upf
+
+NUMERICAL_ORBITAL
+H_gga_6au_60Ry_2s1p.orb
+
+LATTICE_CONSTANT
+1.8897261254578281
+
+LATTICE_VECTORS
+0.000000 1.500000 1.500000
+1.500000 0.000000 1.500000
+1.500000 1.500000 0.000000
+
+ATOMIC_POSITIONS
+Direct
+
+H
+0.0
+2
+0.000000 0.000000 0.000000 0 0 0
+0.500000 0.500000 0.500000 0 0 0
diff --git a/tests/08_RI/63_librpa_soc_symmetry/librpa_producer_manifest.json b/tests/08_RI/63_librpa_soc_symmetry/librpa_producer_manifest.json
new file mode 100644
index 00000000000..55d78ecabd6
--- /dev/null
+++ b/tests/08_RI/63_librpa_soc_symmetry/librpa_producer_manifest.json
@@ -0,0 +1,27 @@
+{
+  "output_dir": "OUT.librpa",
+  "required": [
+    {
+      "pattern": "stru_out.txt",
+      "kind": "stru",
+      "symmetry_rows": 2,
+      "spin_symmetry_rows": 2,
+      "abs_tol": 1e-7
+    },
+    {
+      "pattern": "band_out.txt",
+      "kind": "band"
+    },
+    {
+      "pattern": "bz_sample.txt",
+      "kind": "text"
+    },
+    {
+      "pattern": "KS_wfc_0.dat",
+      "kind": "ks_wfc_v1",
+      "nao_dir": "OUT.autotest",
+      "spinor": true,
+      "abs_tol": 1e-7
+    }
+  ]
+}
diff --git a/tests/integrate/validation_tools/check_librpa_producer.py b/tests/integrate/validation_tools/check_librpa_producer.py
index aebd439e0bd..ee8a47e6e98 100644
--- a/tests/integrate/validation_tools/check_librpa_producer.py
+++ b/tests/integrate/validation_tools/check_librpa_producer.py
@@ -15,6 +15,7 @@
 import struct
 from pathlib import Path
 
+from librpa_structure import check_spin_symmetry
 from librpa_wavefunctions import check_ks_nao, check_velocity, read_ks_wfc
 
 
@@ -94,7 +95,7 @@ def _compare_text(actual, reference, tolerance):
                 )
 
 
-def _check_stru(path, minimum_symmetry_rows):
+def _check_stru(path, minimum_symmetry_rows, spin_symmetry_rows=None):
     lines = path.read_text().splitlines()
     if len(lines) < 8:
         raise ProducerContractError("{} is too short for stru_out".format(path))
@@ -129,6 +130,7 @@ def _check_stru(path, minimum_symmetry_rows):
     if minimum_symmetry_rows is not None:
         if row_match is None or row_match < int(minimum_symmetry_rows):
             raise ProducerContractError("{} has too few symmetry rows".format(path))
+    symmetry_end = None
     if row_match is not None:
         if len(lines) < row_index + 1 + row_match:
             raise ProducerContractError("{} symmetry block is truncated".format(path))
@@ -143,6 +145,12 @@ def _check_stru(path, minimum_symmetry_rows):
                     raise ProducerContractError("{} has a non-integer symmetry rotation".format(path))
             for token in tokens[9:]:
                 _number(token)
+        symmetry_end = row_index + 1 + row_match
+    if spin_symmetry_rows is not None:
+        try:
+            check_spin_symmetry(lines, symmetry_end, path, int(spin_symmetry_rows), row_match)
+        except ValueError as error:
+            raise ProducerContractError(str(error))
 def _check_wfc_nao(path):
     """Validate the text LCAO wavefunction format consumed by GW preprocessing."""
     lines = [line.strip() for line in path.read_text().splitlines() if line.strip()]
@@ -388,7 +396,7 @@ def _check_file(path, entry):
     if kind == "gauge" or kind == "presence":
         return
     if kind == "stru":
-        _check_stru(path, entry.get("symmetry_rows"))
+        _check_stru(path, entry.get("symmetry_rows"), entry.get("spin_symmetry_rows"))
     elif kind == "wfc_nao":
         _check_wfc_nao(path)
     elif kind == "ks_wfc_v1":
@@ -448,7 +456,10 @@ def check_manifest(root, manifest):
                         if dims != (nk, nspin, nbands, nbasis):
                             raise ValueError("{} dimensions disagree with band_out.txt".format(path))
                 if entry.get("kind") == "ks_wfc_v1" and entry.get("nao_dir"):
-                    check_ks_nao(path, root / entry["nao_dir"], float(entry.get("abs_tol", 1.0e-7)))
+                    check_ks_nao(path,
+                                 root / entry["nao_dir"],
+                                 float(entry.get("abs_tol", 1.0e-7)),
+                                 spinor=bool(entry.get("spinor", False)))
             except (ValueError, struct.error) as error:
                 raise ProducerContractError(str(error))
             if entry.get("reference"):
diff --git a/tests/integrate/validation_tools/librpa_structure.py b/tests/integrate/validation_tools/librpa_structure.py
new file mode 100644
index 00000000000..3ab2a0722ca
--- /dev/null
+++ b/tests/integrate/validation_tools/librpa_structure.py
@@ -0,0 +1,36 @@
+"""Validate structure metadata shared by the LibRPA producer checks."""
+
+import math
+
+
+def check_spin_symmetry(lines, symmetry_end, path, minimum_rows, symmetry_rows):
+    if symmetry_end is None or len(lines) <= symmetry_end:
+        raise ValueError("{} has no spin_symmetry block".format(path))
+    header = lines[symmetry_end].split()
+    if len(header) != 3 or header[0] != "spin_symmetry":
+        raise ValueError("{} has an invalid spin_symmetry header".format(path))
+    try:
+        grey_group, spin_source = (int(value) for value in header[1:])
+    except ValueError:
+        raise ValueError("{} has an invalid spin_symmetry header".format(path))
+    if grey_group not in (0, 1) or spin_source not in (0, 1) or minimum_rows < 0:
+        raise ValueError("{} has invalid spin_symmetry flags".format(path))
+    spin_rows = lines[symmetry_end + 1 :]
+    if symmetry_rows is None or symmetry_rows < minimum_rows or len(spin_rows) != symmetry_rows:
+        raise ValueError("{} has an unexpected spin_symmetry row count".format(path))
+    expected_fields = 9 if spin_source == 1 else 1
+    for line in spin_rows:
+        tokens = line.split()
+        if len(tokens) != expected_fields:
+            raise ValueError("{} has an invalid spin_symmetry row".format(path))
+        try:
+            if int(tokens[0]) not in (0, 1):
+                raise ValueError
+        except ValueError:
+            raise ValueError("{} has an invalid spin_symmetry flag".format(path))
+        for token in tokens[1:]:
+            try:
+                if not math.isfinite(float(token)):
+                    raise ValueError
+            except ValueError:
+                raise ValueError("{} has a non-finite spin_symmetry value".format(path))
diff --git a/tests/integrate/validation_tools/librpa_symmetry.py b/tests/integrate/validation_tools/librpa_symmetry.py
new file mode 100644
index 00000000000..8a13646c644
--- /dev/null
+++ b/tests/integrate/validation_tools/librpa_symmetry.py
@@ -0,0 +1,41 @@
+"""Small matrix checks for LibRPA's row-vector symmetry metadata."""
+
+
+def reciprocal_rotation(direct_rotation):
+    """Return the inverse transpose for a direct-space row-vector operation."""
+    if len(direct_rotation) != 9:
+        raise ValueError("Symmetry rotation must have nine elements")
+    a, b, c, d, e, f, g, h, i = map(int, direct_rotation)
+    determinant = a * (e * i - f * h) - b * (d * i - f * g) + c * (d * h - e * g)
+    if abs(determinant) != 1:
+        raise ValueError("Symmetry rotation must be unimodular")
+    return [
+        determinant * (e * i - f * h),
+        determinant * (f * g - d * i),
+        determinant * (d * h - e * g),
+        determinant * (c * h - b * i),
+        determinant * (a * i - c * g),
+        determinant * (b * g - a * h),
+        determinant * (b * f - c * e),
+        determinant * (c * d - a * f),
+        determinant * (a * e - b * d),
+    ]
+
+
+def preserves_bvk_mesh(direct_rotation, mesh, offset=(0.0, 0.0, 0.0), antiunitary=False):
+    """Check that a spatial or antiunitary operation maps an MP mesh to itself."""
+    reciprocal = reciprocal_rotation(direct_rotation)
+    sign = -1 if antiunitary else 1
+    for source in range(3):
+        for target in range(3):
+            numerator = sign * reciprocal[3 * source + target] * mesh[target]
+            if numerator % mesh[source] != 0:
+                return False
+
+    origin = [offset[index] / mesh[index] for index in range(3)]
+    for target in range(3):
+        mapped = sign * sum(origin[source] * reciprocal[3 * source + target] for source in range(3))
+        shift = (mapped - origin[target]) * mesh[target]
+        if abs(shift - round(shift)) > 1e-12:
+            return False
+    return True
diff --git a/tests/integrate/validation_tools/librpa_wavefunctions.py b/tests/integrate/validation_tools/librpa_wavefunctions.py
index b0cdac43f4b..15303a35ef3 100644
--- a/tests/integrate/validation_tools/librpa_wavefunctions.py
+++ b/tests/integrate/validation_tools/librpa_wavefunctions.py
@@ -29,17 +29,22 @@ def read_ks_wfc(path):
     return (nk, nspin, nbands, nbasis), records
 
 
-def check_ks_nao(path, nao_dir, tolerance):
+def check_ks_nao(path, nao_dir, tolerance, spinor=False):
     """Compare the MPI-IO assembly to the independent same-run text writer."""
     (nk, nspin, nbands, nbasis), records = read_ks_wfc(path)
+    if spinor and (nspin != 1 or nbasis % 2):
+        raise ValueError("{} has invalid spinor dimensions".format(path))
     for ik, record in enumerate(records, 1):
         for spin in range(1, nspin + 1):
-            name = "wfk{}_nao.txt".format(ik) if nspin == 1 else "wfk{}s{}_nao.txt".format(ik, spin)
+            if spinor:
+                name = "wfk{}s4_nao.txt".format(ik)
+            else:
+                name = "wfk{}_nao.txt".format(ik) if nspin == 1 else "wfk{}s{}_nao.txt".format(ik, spin)
             source = nao_dir / name
             lines = [line.strip() for line in source.read_text().splitlines() if line.strip()]
             if len(lines) < 4 or "(index of k points)" not in lines[0]:
                 raise ValueError("{} has an invalid text KS header".format(source))
-            text_index = ik + (spin - 1) * nk
+            text_index = ik if spinor else ik + (spin - 1) * nk
             if (int(lines[0].split()[0]), int(lines[2].split()[0]), int(lines[3].split()[0])) != (text_index, nbands, nbasis):
                 raise ValueError("{} dimensions disagree with {}".format(path, source))
             bands = [i for i, line in enumerate(lines) if line.endswith("(band)")]
@@ -51,8 +56,20 @@ def check_ks_nao(path, nao_dir, tolerance):
                 if len(tokens) != 2 * nbasis:
                     raise ValueError("{} has the wrong text coefficient count".format(source))
                 values = [complex(tokens[i], tokens[i + 1]) for i in range(0, len(tokens), 2)]
-                offset = ((spin - 1) * nbands + band) * nbasis
-                for actual, expected in zip(record[offset:offset + nbasis], values):
+                if spinor:
+                    spatial_basis = nbasis // 2
+                    component_values = nbands * spatial_basis
+                    offset = band * spatial_basis
+                    actual_values = [
+                        value
+                        for iw in range(spatial_basis)
+                        for value in (record[offset + iw],
+                                      record[component_values + offset + iw])
+                    ]
+                else:
+                    offset = ((spin - 1) * nbands + band) * nbasis
+                    actual_values = record[offset:offset + nbasis]
+                for actual, expected in zip(actual_values, values):
                     if not math.isfinite(expected.real) or not math.isfinite(expected.imag) or abs(actual - expected) > tolerance:
                         raise ValueError("{} KS payload disagrees with {}".format(path, source))
 
diff --git a/tests/integrate/validation_tools/run_librpa_sym.py b/tests/integrate/validation_tools/run_librpa_sym.py
index 2d855f6d573..61c4606075c 100644
--- a/tests/integrate/validation_tools/run_librpa_sym.py
+++ b/tests/integrate/validation_tools/run_librpa_sym.py
@@ -13,10 +13,13 @@
 import tempfile
 
 from librpa_wavefunctions import check_ks_nao, read_ks_wfc
+from librpa_symmetry import preserves_bvk_mesh
 from run_librpa_mpi import ignore_outputs, mpi_command, use_mpi_solver
 
 
-def check_output(case, mesh, symmetry, nspin):
+def check_output(case, mesh, symmetry, nspin, wfc_nspin=None,
+                 require_spin_symmetry=False, require_antiunitary=False, spinor=False,
+                 offset=(0.0, 0.0, 0.0)):
     out = case / "OUT.librpa"
     rows = [line.split() for line in (out / "bz_sample.txt").read_text().splitlines()]
     assert tuple(map(int, rows[0])) == mesh
@@ -49,13 +52,15 @@ def check_output(case, mesh, symmetry, nspin):
         assert nk == full
     if symmetry == 1 or (symmetry == 0 and mesh == (1, 1, 3)):
         assert nk < full, "The test must exercise an actual reduction"
+    if wfc_nspin is None:
+        wfc_nspin = nspin
     dims, _ = read_ks_wfc(out / "KS_wfc_0.dat")
     bands = (out / "band_out.txt").read_text().split()
-    assert dims[:2] == (nk, nspin)
+    assert dims[:2] == (nk, wfc_nspin)
     assert tuple(map(int, bands[:4])) == (nk, nspin, dims[2], dims[3])
     manifest = json.loads((case / "librpa_producer_manifest.json").read_text())
     tolerance = next(entry["abs_tol"] for entry in manifest["required"] if entry["kind"] == "ks_wfc_v1")
-    check_ks_nao(out / "KS_wfc_0.dat", case / "OUT.autotest", tolerance)
+    check_ks_nao(out / "KS_wfc_0.dat", case / "OUT.autotest", tolerance, spinor=spinor)
     for iq in range(1, nq + 1):
         for rank in (0, 1):
             data = (out / "V_full_{}_r{}.dat".format(iq, rank)).read_bytes()
@@ -67,16 +72,37 @@ def check_output(case, mesh, symmetry, nspin):
         begin = 7 + int(stru[6])
         count, convention = stru[begin].split()
         assert int(count) > 1 and convention == "row"
-        assert len(stru[begin + 1:]) == int(count)
-        assert all(len(row.split()) == 12 for row in stru[begin + 1:])
-        for row in stru[begin + 1:]:
+        row_end = begin + 1 + int(count)
+        assert len(stru) >= row_end
+        assert all(len(row.split()) == 12 for row in stru[begin + 1:row_end])
+        spin_start = row_end
+        antiunitary = [False] * int(count)
+        if require_spin_symmetry or spin_start < len(stru):
+            assert spin_start < len(stru), "SOC stru_out has no spin_symmetry block"
+            header = stru[spin_start].split()
+            assert header[:1] == ["spin_symmetry"] and len(header) == 3
+            assert int(header[1]) in (0, 1)
+            spin_source = int(header[2])
+            assert spin_source in (0, 1)
+            spin_rows = stru[spin_start + 1:]
+            assert len(spin_rows) == int(count)
+            expected_fields = 9 if spin_source == 1 else 1
+            antiunitary = []
+            for row in spin_rows:
+                fields = row.split()
+                assert len(fields) == expected_fields and int(fields[0]) in (0, 1)
+                antiunitary.append(int(fields[0]) == 1)
+                for value in fields[1:]:
+                    assert math.isfinite(float(value))
+            if require_antiunitary:
+                assert any(antiunitary), "shifted magnetic case has no antiunitary operation"
+        for row_index, row in enumerate(stru[begin + 1:row_end]):
             rotation = list(map(int, row.split()[:9]))
-            for i in range(3):
-                for j in range(3):
-                    step = rotation[3 * i + j] * mesh[i] / mesh[j]
-                    assert abs(step - round(step)) < 1e-12, "Operation does not preserve the BvK mesh"
+            assert preserves_bvk_mesh(rotation, mesh, offset, antiunitary[row_index]), (
+                "Operation {} does not preserve the BvK mesh: {}".format(row_index, rotation)
+            )
     return {"case": str(case), "mesh": mesh, "symmetry": symmetry,
-            "nspin": nspin, "nk_full": full, "nk_scf": nk, "nq": nq}
+            "nspin": nspin, "nspin_wfc": wfc_nspin, "nk_full": full, "nk_scf": nk, "nq": nq}
 
 
 def run(args, root):
@@ -93,7 +119,7 @@ def run(args, root):
                ABACUS_MPIEXEC_PREFLAGS=" ".join(args.mpi_preflag),
                ABACUS_MPIEXEC_POSTFLAGS=" ".join(args.mpi_postflag))
     cases = (source / "CASES_LIBRPA_PRODUCER.txt").read_text().split()
-    for name, nspin in zip(cases, (1, 2)):
+    for name, nspin in zip(cases[:2], (1, 2)):
         for mesh in ((1, 1, 3), (2, 2, 2)):
             for symmetry in (-1, 0, 1):
                 case = work / "{}_{}_s{}".format(name, "".join(map(str, mesh)), symmetry)
@@ -115,8 +141,63 @@ def run(args, root):
                 result = check_output(case, mesh, symmetry, nspin)
                 results.append(result)
                 print(json.dumps(result), flush=True)
+    soc_case = work / "63_librpa_soc_symmetry_222_s1"
+    shutil.copytree(source / "63_librpa_soc_symmetry", soc_case, ignore=ignore_outputs)
+    use_mpi_solver(soc_case)
+    soc_input = soc_case / "INPUT"
+    soc_input.write_text(re.sub(r"(?m)^symmetry\s+.*$", "symmetry 1", soc_input.read_text()))
+    with soc_input.open("a") as stream:
+        stream.write("scf_thr 1e-10\n")
+    (soc_case / "KPT").write_text("K_POINTS\n0\nGamma\n2 2 2 0 0 0\n")
+    with (soc_case / "producer.log").open("w") as log:
+        subprocess.run(mpi_command(args, 2), cwd=soc_case, env=env,
+                       stdout=log, stderr=subprocess.STDOUT, check=True, timeout=180)
+    soc_result = check_output(soc_case, (2, 2, 2), 1, 1,
+                              require_spin_symmetry=True, spinor=True)
+    results.append(soc_result)
+    print(json.dumps(soc_result), flush=True)
+
+    magnetic_case = work / "63_librpa_soc_symmetry_shifted_magnetic"
+    shutil.copytree(source / "63_librpa_soc_symmetry", magnetic_case, ignore=ignore_outputs)
+    use_mpi_solver(magnetic_case)
+    magnetic_stru = []
+    diagonal_vectors = iter(("3.000000 0.000000 0.000000",
+                             "0.000000 3.000000 0.000000",
+                             "0.000000 0.000000 3.000000"))
+    replace_vectors = False
+    vector_index = 0
+    for line in (magnetic_case / "STRU").read_text().splitlines():
+        if line == "LATTICE_VECTORS":
+            replace_vectors = True
+            magnetic_stru.append(line)
+            continue
+        if replace_vectors:
+            line = next(diagonal_vectors)
+            vector_index += 1
+            if vector_index == 3:
+                replace_vectors = False
+        fields = line.split()
+        if len(fields) == 6 and fields[3:] == ["0", "0", "0"]:
+            line += " mag 0 0 1"
+        magnetic_stru.append(line)
+    (magnetic_case / "STRU").write_text("\n".join(magnetic_stru) + "\n")
+    magnetic_input = magnetic_case / "INPUT"
+    magnetic_input.write_text(re.sub(r"(?m)^symmetry\s+.*$", "symmetry 1", magnetic_input.read_text()))
+    with magnetic_input.open("a") as stream:
+        stream.write("scf_thr 1e-10\n")
+    (magnetic_case / "KPT").write_text("K_POINTS\n0\nGamma\n2 3 3 0.25 0 0\n")
+    with (magnetic_case / "producer.log").open("w") as log:
+        subprocess.run(mpi_command(args, 2), cwd=magnetic_case, env=env,
+                       stdout=log, stderr=subprocess.STDOUT, check=True, timeout=180)
+    magnetic_result = check_output(magnetic_case, (2, 3, 3), 1, 1,
+                                   require_spin_symmetry=True,
+                                   require_antiunitary=True,
+                                   spinor=True,
+                                   offset=(0.25, 0.0, 0.0))
+    results.append(magnetic_result)
+    print(json.dumps(magnetic_result), flush=True)
     (root / "summary.json").write_text(json.dumps(results, indent=2) + "\n")
-    print("Symmetry producer checks: PASS (12 two-rank runs)")
+    print("Symmetry producer checks: PASS (12 baseline runs and two SOC runs)")
 
 
 def main():
diff --git a/tests/integrate/validation_tools/run_librpa_wavefunction_checks.sh b/tests/integrate/validation_tools/run_librpa_wavefunction_checks.sh
index ed0d6316f72..876851f16e0 100755
--- a/tests/integrate/validation_tools/run_librpa_wavefunction_checks.sh
+++ b/tests/integrate/validation_tools/run_librpa_wavefunction_checks.sh
@@ -7,4 +7,4 @@ if [ "$#" -ne 1 ]; then
     exit 2
 fi
 
-python3 -m unittest discover -s "$1" -p test_librpa_wavefunctions.py
+python3 -m unittest discover -s "$1" -p 'test_librpa_*.py'
diff --git a/tests/integrate/validation_tools/test_librpa_symmetry.py b/tests/integrate/validation_tools/test_librpa_symmetry.py
new file mode 100644
index 00000000000..6377e3e9a2a
--- /dev/null
+++ b/tests/integrate/validation_tools/test_librpa_symmetry.py
@@ -0,0 +1,27 @@
+"""Check reciprocal-grid symmetry handling for reader-v1 metadata."""
+
+import unittest
+
+from librpa_symmetry import preserves_bvk_mesh, reciprocal_rotation
+
+
+class LibrpaSymmetryTests(unittest.TestCase):
+    def test_row_rotation_uses_inverse_transpose_for_k_points(self):
+        direct = [1, 0, 0, 1, -1, 0, 1, 0, -1]
+        self.assertEqual(reciprocal_rotation(direct), [1, 1, 1, 0, -1, 0, 0, 0, -1])
+        self.assertTrue(preserves_bvk_mesh(direct, (1, 1, 3)))
+        self.assertFalse(preserves_bvk_mesh(direct, (2, 1, 3)))
+
+    def test_antiunitary_action_preserves_shifted_mesh(self):
+        inversion = [-1, 0, 0, 0, -1, 0, 0, 0, -1]
+        offset = (0.25, 0.0, 0.0)
+        self.assertFalse(preserves_bvk_mesh(inversion, (2, 3, 3), offset))
+        self.assertTrue(preserves_bvk_mesh(inversion, (2, 3, 3), offset, antiunitary=True))
+
+    def test_rejects_non_unimodular_rotation(self):
+        with self.assertRaisesRegex(ValueError, "unimodular"):
+            reciprocal_rotation([2, 0, 0, 0, 1, 0, 0, 0, 1])
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tests/integrate/validation_tools/test_librpa_wavefunctions.py b/tests/integrate/validation_tools/test_librpa_wavefunctions.py
index 0344400d53b..51a9a3c4d60 100644
--- a/tests/integrate/validation_tools/test_librpa_wavefunctions.py
+++ b/tests/integrate/validation_tools/test_librpa_wavefunctions.py
@@ -55,6 +55,67 @@ def test_spin_blocks_and_global_text_indices(self):
         with self.assertRaisesRegex(ValueError, "disagrees"):
             check_ks_nao(self.ks, self.root, 1e-7)
 
+    def test_spinor_text_file_uses_spinor_tag(self):
+        # Binary v1 stores component-major blocks; the text writer interleaves
+        # the two spinor components at each spatial orbital.
+        values = (1, 0.5, 2, -0.5, 3, 0.25, 4, -0.25,
+                  5, 0.125, 6, -0.125, 7, 0.75, 8, -0.75)
+        header = struct.pack("<6iiq", -12345679, 28, 1, 1, 2, 4, 1, 36)
+        payload = struct.pack("<16d", *values)
+        self.ks.write_bytes(header + payload)
+        (self.root / "wfk1s4_nao.txt").write_text(
+            "1 (index of k points)\n0 0 0\n2 (number of bands)\n4 (number of orbitals)\n"
+            "1 (band)\n0 (Ry)\n1 (Occupations)\n"
+            "1 0.5 5 0.125 2 -0.5 6 -0.125\n"
+            "2 (band)\n1 (Ry)\n0 (Occupations)\n"
+            "3 0.25 7 0.75 4 -0.25 8 -0.75\n")
+        check_ks_nao(self.ks, self.root, 1e-7, spinor=True)
+
+    def test_manifest_propagates_soc_contract_requirements(self):
+        output = self.root / "OUT.librpa"
+        output.mkdir()
+        values = (1, 0.5, 2, -0.5, 3, 0.25, 4, -0.25,
+                  5, 0.125, 6, -0.125, 7, 0.75, 8, -0.75)
+        output.joinpath("KS_wfc_0.dat").write_bytes(
+            struct.pack("<6iiq", -12345679, 28, 1, 1, 2, 4, 1, 36)
+            + struct.pack("<16d", *values)
+        )
+        output.joinpath("band_out.txt").write_text("1\n1\n2\n4\n0\n")
+        output.joinpath("stru_out.txt").write_text(
+            "1 0 0\n0 1 0\n0 0 1\n1 0 0\n0 1 0\n0 0 1\n"
+            "1\n0 0 0 0\n2 row\n"
+            "1 0 0 0 1 0 0 0 1 0 0 0\n"
+            "1 0 0 0 1 0 0 0 1 0 0 0\n"
+            "spin_symmetry 0 1\n"
+            "0 1 0 0 0 0 0 0 1\n"
+            "1 1 0 0 0 0 0 0 1\n"
+        )
+        nao = self.root / "OUT.autotest"
+        nao.mkdir()
+        nao.joinpath("wfk1s4_nao.txt").write_text(
+            "1 (index of k points)\n0 0 0\n2 (number of bands)\n4 (number of orbitals)\n"
+            "1 (band)\n0 (Ry)\n1 (Occupations)\n"
+            "1 0.5 5 0.125 2 -0.5 6 -0.125\n"
+            "2 (band)\n1 (Ry)\n0 (Occupations)\n"
+            "3 0.25 7 0.75 4 -0.25 8 -0.75\n"
+        )
+        manifest = {
+            "output_dir": "OUT.librpa",
+            "required": [
+                {"pattern": "stru_out.txt", "kind": "stru", "symmetry_rows": 2,
+                 "spin_symmetry_rows": 2},
+                {"pattern": "band_out.txt", "kind": "band"},
+                {"pattern": "KS_wfc_0.dat", "kind": "ks_wfc_v1", "nao_dir": "OUT.autotest",
+                 "spinor": True},
+            ],
+        }
+        check_manifest(self.root, manifest)
+        output.joinpath("stru_out.txt").write_text(
+            output.joinpath("stru_out.txt").read_text().split("spin_symmetry", 1)[0]
+        )
+        with self.assertRaisesRegex(ProducerContractError, "spin_symmetry"):
+            check_manifest(self.root, manifest)
+
     def test_truncated_header_and_payload(self):
         for data in (self.data[:20], self.data[:-1], self.data + b"x"):
             self.ks.write_bytes(data)