diff --git a/.github/workflows/build_test_makefile.yml b/.github/workflows/build_test_makefile.yml
index bb6f25c206d..501f3247531 100644
--- a/.github/workflows/build_test_makefile.yml
+++ b/.github/workflows/build_test_makefile.yml
@@ -17,17 +17,34 @@ jobs:
           build_args: "CXX=mpiicpx ELPA_LIB_DIR=/usr/local/lib ELPA_INCLUDE_DIR=/usr/local/include CEREAL_DIR=/usr/include/cereal OPENMP=ON"
           name: "Build with Makefile & Intel compilers"
     name: ${{ matrix.name }}
-    container: ghcr.io/deepmodeling/abacus-${{ matrix.tag }}
     steps:
       - name: Checkout
         uses: actions/checkout@v7.0.1
+      # A job-level container is pulled before steps run, so cleanup must happen
+      # on the hosted runner before starting the large Intel toolchain image.
+      - name: Free disk space before pulling the build image
+        run: |
+          df -h /
+          sudo rm -rf /usr/share/dotnet /usr/local/lib/android /opt/ghc
+          docker image prune --all --force
+          df -h /
       - name: Build
-        # TODO: Remove nlohmann-json installation here since Dockerfile will have installed it.
+        env:
+          BUILD_IMAGE: ghcr.io/deepmodeling/abacus-${{ matrix.tag }}
+          BUILD_ARGS: ${{ matrix.build_args }}
+        # Keep the same image and build command as the container job.
         run: |
-          apt update && apt install -y nlohmann-json3-dev
-          export I_MPI_CXX=icpx
-          chmod a+x generate_build_info.sh
-          cd source
-          mkdir build && cd build
-          make -f ../Makefile -j2 ${{ matrix.build_args }}
-
+          docker run --rm \
+            --volume "$GITHUB_WORKSPACE:/workspace" \
+            --workdir /workspace \
+            --env BUILD_ARGS \
+            --env I_MPI_CXX=icpx \
+            --entrypoint /bin/bash \
+            "$BUILD_IMAGE" -e -c '
+              git config --global --add safe.directory /workspace
+              apt update && apt install -y nlohmann-json3-dev
+              chmod a+x generate_build_info.sh
+              cd source
+              mkdir -p build && cd build
+              make -f ../Makefile -j2 $BUILD_ARGS
+            '
diff --git a/source/source_base/module_out/csr_reader.cpp b/source/source_base/module_out/csr_reader.cpp
index 7697f057bbf..0345c0723c0 100644
--- a/source/source_base/module_out/csr_reader.cpp
+++ b/source/source_base/module_out/csr_reader.cpp
@@ -1,4 +1,7 @@
 #include "csr_reader.h"
+#include <cmath>
+#include <set>
+#include <limits>
 #include "source_base/tool_quit.h"
 
 namespace ModuleIO
@@ -11,180 +14,187 @@ csrFileReader<T>::csrFileReader(const std::string& filename) : FileReader(filena
     parseFile();
 }
 
-// function to parse file
+// Format-specific metadata is followed by the same CSR block reader.
 template <typename T>
 void csrFileReader<T>::parseFile()
 {
-    // Check if file is open
-    if (!isOpen())
+    RCoordinates.clear();
+    sparse_matrices.clear();
+    nextDataLine();
+    const std::string line = ss.str();
+    const std::size_t begin = line.find_first_not_of(" \t");
+    if (line.compare(begin, 5, "STEP:") == 0)
     {
-        ModuleBase::WARNING_QUIT("csrFileReader::parseFile", "File is not open");
+        parseCompactHeader();
     }
+    else
+    {
+        parseCellHeader();
+    }
+    if (matrixDimension <= 0 || numberOfR < 0)
+    {
+        ModuleBase::WARNING_QUIT("csrFileReader::parseFile", "Invalid CSR dimensions or R-block count");
+    }
+    readBlocks();
+    validateFrameEnd();
+}
 
-    std::string tmp_string;
-
-    // Read the step
-    readLine();
-    ss >> tmp_string >> tmp_string >> tmp_string >> step;
-
-    //std::cout << " step is " << step << std::endl; 
-    // Read the title
-    readLine();
-    // Read the total spin
-    readLine();
-    // Read the spin index
-    readLine();
+template <typename T>
+void csrFileReader<T>::parseCompactHeader()
+{
+    const auto invalid = []() {
+        ModuleBase::WARNING_QUIT("csrFileReader::parseCompactHeader", "Malformed compact CSR file");
+    };
+    std::string label;
+    if (!(ss >> label >> step) || label != "STEP:") { invalid(); }
+    const auto header_integer = [&](const std::string& prefix, int& value) {
+        readLine();
+        const std::string line = ss.str();
+        const std::size_t begin = line.find_first_not_of(" \t");
+        const std::size_t colon = line.find(':');
+        if (begin == std::string::npos || line.compare(begin, prefix.size(), prefix) != 0
+            || colon == std::string::npos) { invalid(); }
+        std::istringstream values(line.substr(colon + 1));
+        if (!(values >> value)) { invalid(); }
+    };
+    header_integer("Matrix Dimension of ", matrixDimension);
+    header_integer("Matrix number of ", numberOfR);
+}
 
-    // Read the matrix dimension
+// Cell-aware headers used by HContainer writers include spin and geometry.
+template <typename T>
+void csrFileReader<T>::parseCellHeader()
+{
+    const auto invalid = []() {
+        ModuleBase::WARNING_QUIT("csrFileReader::parseCellHeader", "Malformed cell-aware CSR header");
+    };
+    std::string marker;
+    std::string ionic;
+    std::string step_label;
+    if (!(ss >> marker >> ionic >> step_label >> step)
+        || marker != "---" || ionic != "Ionic" || step_label != "Step") { invalid(); }
+    readLine(); // matrix description
     readLine();
-    ss >> matrixDimension;
-    // std::cout << " mat dim is " << matrixDimension << std::endl; 
-
-    // Read the number of R
+    int spins = 0;
+    if (!(ss >> spins) || spins <= 0) { invalid(); }
     readLine();
-    ss >> numberOfR;
-
-    // Read cell
-    // Skip empty line before ucell info if present
+    int spin = 0;
+    if (!(ss >> spin) || spin < 1 || spin > spins) { invalid(); }
     readLine();
-    std::string line = ss.str();
-    if (line.empty() || line.find_first_not_of(" \t\n\r") == std::string::npos)
-    {
-        // Empty line, read next line for latName
-        readLine();
-    }
-    // Now ss contains latName (or the line after empty line)
-    // We don't need to use latName, just continue reading
-    
-    // Read lat0, latvec (3 lines), atom labels (6 lines total including latName)
-    readLine(); // lat0
-    readLine(); // latvec e1
-    readLine(); // latvec e2
-    readLine(); // latvec e3
-    readLine(); // atom labels
-    
-    // Read atom numbers
+    if (!(ss >> matrixDimension)) { invalid(); }
     readLine();
-    std::istringstream iss(ss.str());
-    int natom = 0;
-    int total_natom = 0;
-    while (iss >> natom)
+    if (!(ss >> numberOfR)) { invalid(); }
+    nextDataLine(); // lattice name, possibly after an empty line
+    for (int i = 0; i < 5; ++i) { readLine(); } // scale, three vectors, species
+    readLine(); // atom counts
+    int atoms = 0;
+    int count = 0;
+    while (true)
     {
-        total_natom += natom;
-    }
-    
-    // Read "Direct" line and atom coordinates
-    for (int i = 0; i < total_natom + 1; i++)
-    {
-        readLine(); // Direct + atom coordinates
+        ss >> std::ws;
+        if (ss.eof() || ss.peek() == '#') { break; }
+        if (!(ss >> count) || count < 0 || count > std::numeric_limits<int>::max() - atoms)
+        { invalid(); }
+        atoms += count;
     }
+    if (atoms <= 0) { invalid(); }
+    readLine(); // coordinate convention
+    for (int i = 0; i < atoms; ++i) { readLine(); }
+}
 
-    // Skip empty lines and CSR format comment block
-    // (lines starting with '#' or containing only whitespace)
-    // Use readLine() since ifs is private in FileReader
-    bool found_data = false;
-    while (!found_data)
+// Skip only whitespace and comments, never malformed numerical data.
+template <typename T>
+void csrFileReader<T>::nextDataLine()
+{
+    do
     {
         readLine();
-        std::string line = ss.str();
-        size_t start = line.find_first_not_of(" \t");
-        if (start != std::string::npos && line[start] != '#')
-        {
-            // This is a data line (first R-block), ss already holds it
-            found_data = true;
-        }
-    }
-
-    // Read the matrices
-    // Note: ss already contains the first R-block line from the skip loop above
-    for (int i = 0; i < numberOfR; i++)
-    {
-        // std::cout << " read R " << i+1 << std::endl;
+        ss >> std::ws;
+    } while (ss.eof() || ss.peek() == '#');
+}
 
-        std::vector<int> RCoord(3);
-        int nonZero = 0;
+template <typename T>
+void csrFileReader<T>::nextDataToken()
+{
+    ss >> std::ws;
+    if (ss.eof() || ss.peek() == '#') { nextDataLine(); }
+}
 
-        if (i > 0)
+template <typename T>
+void csrFileReader<T>::readBlocks()
+{
+    const auto invalid = []() {
+        ModuleBase::WARNING_QUIT("csrFileReader::readBlocks", "Malformed CSR block");
+    };
+    std::set<std::vector<int>> seen;
+    for (int ir = 0; ir < numberOfR; ++ir)
+    {
+        nextDataLine();
+        std::vector<int> coordinate(3);
+        int nnz = 0;
+        if (!(ss >> coordinate[0] >> coordinate[1] >> coordinate[2] >> nnz) || nnz < 0)
+        { invalid(); }
+        const std::size_t dimension = matrixDimension;
+        if (static_cast<std::size_t>(nnz) > dimension * dimension || !seen.insert(coordinate).second)
+        { invalid(); }
+        std::vector<T> values(nnz);
+        std::vector<int> columns(nnz);
+        std::vector<int> pointers(dimension + 1);
+        for (int i = 0; i < nnz; ++i)
         {
-            // For subsequent R-blocks, skip empty lines to find next R-coordinate line
-            bool found = false;
-            while (!found)
-            {
-                readLine();
-                std::string line = ss.str();
-                size_t start = line.find_first_not_of(" \t");
-                if (start != std::string::npos && line[start] != '#')
-                {
-                    found = true;
-                }
-            }
+            nextDataToken();
+            if (!(ss >> values[i]) || !std::isfinite(std::abs(values[i]))) { invalid(); }
         }
-        ss >> RCoord[0] >> RCoord[1] >> RCoord[2] >> nonZero;
-        RCoordinates.push_back(RCoord);
-
-        std::vector<T> csr_values(nonZero);
-        std::vector<int> csr_col_ind(nonZero);
-        std::vector<int> csr_row_ptr(matrixDimension + 1);
-
-        // read CSR values
-        readLine();
-        // std::cout << " ss1: " << ss.str() << std::endl;
-
-        readLine();
-    size_t count1 = 0;
-        while (count1 < nonZero)
+        for (int i = 0; i < nnz; ++i)
         {
-            if (ss.eof() || ss.fail())
-            {
-                readLine();
+            nextDataToken();
+            if (!(ss >> columns[i]) || columns[i] < 0 || columns[i] >= matrixDimension) { invalid(); }
         }
-            if (ss >> csr_values[count1])
-            {
-                count1++;
-            }
-    }
-        // std::cout << "count1=" << count1 << std::endl;
-
-        // read CSR column indices
-        readLine();
-        // std::cout << " ss2: " << ss.str() << std::endl;
-
-    size_t count2 = 0;
-        while (count2 < nonZero)
+        for (int i = 0; i <= matrixDimension; ++i)
         {
-            if (ss.eof() || ss.fail())
-            {
-                readLine();
+            nextDataToken();
+            if (!(ss >> pointers[i]) || pointers[i] < 0 || pointers[i] > nnz) { invalid(); }
+            if (i > 0 && pointers[i] < pointers[i - 1]) { invalid(); }
         }
-            if (ss >> csr_col_ind[count2])
+        if (pointers.front() != 0 || pointers.back() != nnz) { invalid(); }
+        ss >> std::ws;
+        if (!ss.eof() && ss.peek() != '#') { invalid(); }
+        for (int row = 0; row < matrixDimension; ++row)
+        {
+            for (int i = pointers[row] + 1; i < pointers[row + 1]; ++i)
             {
-                count2++;
+                if (columns[i] <= columns[i - 1]) { invalid(); }
             }
+        }
+        SparseMatrix<T> matrix(matrixDimension, matrixDimension);
+        matrix.readCSR(values, columns, pointers);
+        RCoordinates.push_back(coordinate);
+        sparse_matrices.push_back(matrix);
     }
-        // std::cout << "count2=" << count2 << std::endl;
-
-        // read row pointers
-        readLine();
-        // std::cout << " ss3: " << ss.str() << std::endl;
+}
 
-    size_t count3 = 0;
-        while (count3 < matrixDimension + 1)
+// Only blank/comment lines, EOF or an unconsumed next frame may follow a frame.
+template <typename T>
+void csrFileReader<T>::validateFrameEnd()
+{
+    std::string line;
+    while (peekLine(line))
+    {
+        std::istringstream tokens(line);
+        tokens >> std::ws;
+        if (tokens.eof() || tokens.peek() == '#')
         {
-            if (ss.eof() || ss.fail())
-            {
-                readLine();
+            readLine();
+            continue;
         }
-            if (ss >> csr_row_ptr[count3])
-            {
-                count3++;
-            }
-    }
-        // std::cout << "count3=" << count3 << std::endl;
-
-        // create sparse matrix
-        SparseMatrix<T> matrix(matrixDimension, matrixDimension);
-        matrix.readCSR(csr_values, csr_col_ind, csr_row_ptr);
-        sparse_matrices.push_back(matrix);
+        std::string marker;
+        tokens >> marker;
+        if (marker == "STEP:") { return; }
+        std::string ionic;
+        std::string step_label;
+        if (marker == "---" && tokens >> ionic >> step_label
+            && ionic == "Ionic" && step_label == "Step") { return; }
+        ModuleBase::WARNING_QUIT("csrFileReader::validateFrameEnd", "Unexpected data after CSR frame");
     }
 }
 
diff --git a/source/source_base/module_out/csr_reader.h b/source/source_base/module_out/csr_reader.h
index 6b38db07f3d..8214fc33981 100644
--- a/source/source_base/module_out/csr_reader.h
+++ b/source/source_base/module_out/csr_reader.h
@@ -12,7 +12,12 @@ namespace ModuleIO
 /**
  * @brief Class to read CSR file
  * @details This class is used to read CSR file
- * the default format is like:
+ * Both save_sparse's compact three-line header (STEP, Matrix Dimension,
+ * Matrix number) and the cell-aware format below are accepted. Compact CSR
+ * blocks contain an R coordinate/nnz line, values, zero-based column indices,
+ * and dimension+1 row pointers. Invalid dimensions, indices and truncated data
+ * are rejected for both formats rather than producing an empty matrix.
+ * The cell-aware format is like:
  * ```
  *  --- Ionic Step 1 ---
  *  # print density matrix in real space DM(R)
@@ -72,11 +77,17 @@ class csrFileReader : public FileReader
     int getMatrixDimension() const;
 
   private:
+    void parseCompactHeader();
+    void parseCellHeader();
+    void readBlocks();
+    void validateFrameEnd();
+    void nextDataLine();
+    void nextDataToken();
     std::vector<std::vector<int>> RCoordinates;
     std::vector<SparseMatrix<T>> sparse_matrices;
-    int step;
-    int matrixDimension;
-    int numberOfR;
+    int step = 0;
+    int matrixDimension = 0;
+    int numberOfR = 0;
 };
 
 } // namespace ModuleIO
diff --git a/source/source_base/module_out/file_reader.cpp b/source/source_base/module_out/file_reader.cpp
index 7211ac44456..f8f5ee35ced 100644
--- a/source/source_base/module_out/file_reader.cpp
+++ b/source/source_base/module_out/file_reader.cpp
@@ -47,6 +47,16 @@ void FileReader::readLine()
     }
 }
 
+bool FileReader::peekLine(std::string& line)
+{
+    const std::streampos position = ifs.tellg();
+    if (position == std::streampos(-1)) { return false; }
+    const bool found = static_cast<bool>(std::getline(ifs, line));
+    ifs.clear();
+    ifs.seekg(position);
+    return found;
+}
+
 void FileReader::read_ucell()
 {
     if (ifs.eof())
diff --git a/source/source_base/module_out/file_reader.h b/source/source_base/module_out/file_reader.h
index 8f484c39693..4c4dc679a51 100644
--- a/source/source_base/module_out/file_reader.h
+++ b/source/source_base/module_out/file_reader.h
@@ -33,6 +33,10 @@ class FileReader
 
     std::stringstream ss;
 
+  protected:
+    // Inspect the next line without consuming it; false at end of file.
+    bool peekLine(std::string& line);
+
   private:
     std::ifstream ifs;
 };
diff --git a/source/source_base/test/CMakeLists.txt b/source/source_base/test/CMakeLists.txt
index c9c0a3f8500..1dc92e4e4cc 100644
--- a/source/source_base/test/CMakeLists.txt
+++ b/source/source_base/test/CMakeLists.txt
@@ -268,7 +268,7 @@ addtest(
 addtest(
   TARGET MODULE_BASE_csr_reader
   LIBS parameter base device
-  SOURCES csr_reader_test.cpp)
+  SOURCES csr_reader_test.cpp test_csr_reader.cpp)
 
 if(ENABLE_GOOGLEBENCH)
   addtest(
diff --git a/source/source_base/test/test_csr_reader.cpp b/source/source_base/test/test_csr_reader.cpp
new file mode 100644
index 00000000000..7f805b886ab
--- /dev/null
+++ b/source/source_base/test/test_csr_reader.cpp
@@ -0,0 +1,151 @@
+#include "source_base/module_out/csr_reader.h"
+#include "gtest/gtest.h"
+#include <cstdio>
+#include <fstream>
+#include <complex>
+#include <iterator>
+
+namespace
+{
+class CsrFormats : public testing::TestWithParam<bool>
+{
+protected:
+    std::string path;
+    void SetUp() override { path = testing::TempDir() + "/abacus_compact_csr_test.csr"; }
+    void TearDown() override { std::remove(path.c_str()); }
+    void write(const std::string& blocks, const std::string& dimension, const std::string& count)
+    {
+        std::ofstream out(path);
+        if (GetParam())
+        {
+            out << "STEP: 7\nMatrix Dimension of Hexxs_0(R): " << dimension
+                << "\nMatrix number of Hexxs_0(R): " << count << '\n';
+        }
+        else
+        {
+            out << " --- Ionic Step 7 ---\n # matrix\n1 # spins\n1 # spin\n"
+                << dimension << " # basis\n" << count << " # R blocks\n\n"
+                << "user_defined_lattice\n1\n1 0 0\n0 1 0\n0 0 1\nSi\n1\nDirect\n0 0 0\n";
+        }
+        out << blocks;
+    }
+};
+TEST_P(CsrFormats, RealMultipleRAndWrappedValues)
+{
+    write("0 0 0 2\n 1.25\n -2\n 0 1\n 0 1 2\n0 0 1 1\n 3\n 1\n 0 1 1\n", "2", "2");
+    ModuleIO::csrFileReader<double> reader(path);
+    EXPECT_EQ(reader.getStep(), 7);
+    EXPECT_EQ(reader.getMatrixDimension(), 2);
+    EXPECT_EQ(reader.getNumberOfR(), 2);
+    EXPECT_DOUBLE_EQ(reader.getMatrix(0)(0, 0), 1.25);
+    EXPECT_DOUBLE_EQ(reader.getMatrix(0)(1, 1), -2);
+    EXPECT_DOUBLE_EQ(reader.getMatrix(0, 0, 1)(0, 1), 3);
+}
+TEST_P(CsrFormats, ComplexValues)
+{
+    write("0 0 0 2\n (1,2) (-3,-4)\n0 1\n0 1 2\n", "2", "1");
+    ModuleIO::csrFileReader<std::complex<double>> reader(path);
+    EXPECT_EQ(reader.getMatrix(0)(0, 0), std::complex<double>(1, 2));
+    EXPECT_EQ(reader.getMatrix(0)(1, 1), std::complex<double>(-3, -4));
+}
+TEST_P(CsrFormats, EmptyRankMatrix)
+{
+    write("", "2", "0");
+    ModuleIO::csrFileReader<double> reader(path);
+    EXPECT_EQ(reader.getMatrixDimension(), 2);
+    EXPECT_EQ(reader.getNumberOfR(), 0);
+}
+TEST_P(CsrFormats, RejectsInvalidDimension)
+{
+    write("", "0", "0");
+    EXPECT_EXIT({ModuleIO::csrFileReader<double> reader(path);}, testing::ExitedWithCode(1), "");
+}
+TEST_P(CsrFormats, RejectsMalformedHeader)
+{
+    write("", "wrong", "1");
+    EXPECT_EXIT({ModuleIO::csrFileReader<double> reader(path);}, testing::ExitedWithCode(1), "");
+}
+TEST_P(CsrFormats, RejectsTruncatedValues)
+{
+    write("0 0 0 2\n1\n", "2", "1");
+    EXPECT_EXIT({ModuleIO::csrFileReader<double> reader(path);}, testing::ExitedWithCode(1), "");
+}
+TEST_P(CsrFormats, RejectsInvalidColumn)
+{
+    write("0 0 0 1\n1\n2\n0 1 1\n", "2", "1");
+    EXPECT_EXIT({ModuleIO::csrFileReader<double> reader(path);}, testing::ExitedWithCode(1), "");
+}
+TEST_P(CsrFormats, RejectsInvalidRowPointers)
+{
+    write("0 0 0 1\n1\n0\n0 1 0\n", "2", "1");
+    EXPECT_EXIT({ModuleIO::csrFileReader<double> reader(path);}, testing::ExitedWithCode(1), "");
+}
+TEST_P(CsrFormats, RejectsMalformedNumber)
+{
+    write("0 0 0 1\nbad\n0\n0 1 1\n", "2", "1");
+    EXPECT_EXIT({ModuleIO::csrFileReader<double> reader(path);}, testing::ExitedWithCode(1), "");
+}
+
+TEST_P(CsrFormats, CommentsAndEmptyBlock)
+{
+    write("# blocks\n0 0 0 0\n# values\n\n# columns\n\n# pointers\n0 0 0\n"
+          "\n0 0 1 1\n# values\n2.5\n# columns\n1\n# pointers\n0 1 1\n", "2", "2");
+    ModuleIO::csrFileReader<double> reader(path);
+    EXPECT_TRUE(reader.getMatrix(0).getElements().empty());
+    EXPECT_DOUBLE_EQ(reader.getMatrix(1)(0, 1), 2.5);
+}
+TEST_P(CsrFormats, RejectsDuplicateR)
+{
+    write("0 0 0 1\n1\n0\n0 1 1\n0 0 0 1\n2\n0\n0 1 1\n", "2", "2");
+    EXPECT_EXIT({ModuleIO::csrFileReader<double> reader(path);}, testing::ExitedWithCode(1), "");
+}
+TEST_P(CsrFormats, RejectsExtraValueShiftingRowPointers)
+{
+    write("0 0 0 1\n1 0\n0\n0 1 1\n", "2", "1");
+    EXPECT_EXIT({ModuleIO::csrFileReader<double> reader(path);}, testing::ExitedWithCode(1), "");
+}
+TEST_P(CsrFormats, AcceptsCommentAfterFinalRowPointer)
+{
+    write("0 0 0 1\n1\n0\n0 1 1 # end of block\n", "2", "1");
+    ModuleIO::csrFileReader<double> reader(path);
+    EXPECT_DOUBLE_EQ(reader.getMatrix(0)(0, 0), 1.0);
+}
+TEST_P(CsrFormats, RejectsWrappedLeftoverRowPointer)
+{
+    write("0 0 0 1\n1 0\n0\n0 1\n1\n", "2", "1");
+    EXPECT_EXIT({ModuleIO::csrFileReader<double> reader(path);}, testing::ExitedWithCode(1), "");
+}
+TEST_P(CsrFormats, PreservesNextFrameAfterComments)
+{
+    write("0 0 0 1\n1\n0\n0\n1 1\n\n# frame end\n", "2", "1");
+    std::ifstream input(path);
+    const std::istreambuf_iterator<char> begin(input);
+    const std::istreambuf_iterator<char> end;
+    std::string frame(begin, end);
+    input.close();
+    const std::size_t step_position = frame.find('7');
+    ASSERT_NE(step_position, std::string::npos);
+    frame.replace(step_position, 1, "8");
+    const std::string first_block = "0 0 0 1\n1\n";
+    const std::string second_block = "0 0 0 1\n2.5\n";
+    const std::size_t block_position = frame.find(first_block);
+    const std::size_t block_size = first_block.size();
+    ASSERT_NE(block_position, std::string::npos);
+    frame.replace(block_position, block_size, second_block);
+    std::ofstream append(path, std::ios::app);
+    append << frame;
+    append.close();
+    ModuleIO::csrFileReader<double> reader(path);
+    EXPECT_DOUBLE_EQ(reader.getMatrix(0)(0, 0), 1.0);
+    EXPECT_EQ(reader.getStep(), 7);
+    reader.parseFile();
+    EXPECT_EQ(reader.getStep(), 8);
+    EXPECT_DOUBLE_EQ(reader.getMatrix(0)(0, 0), 2.5);
+}
+TEST_P(CsrFormats, RejectsDataAfterTrailingComment)
+{
+    write("0 0 0 1\n1\n0\n0 1 1\n# end\n\n9\n", "2", "1");
+    EXPECT_EXIT({ModuleIO::csrFileReader<double> reader(path);}, testing::ExitedWithCode(1), "");
+}
+INSTANTIATE_TEST_SUITE_P(BothHeaders, CsrFormats, testing::Bool());
+}
diff --git a/source/source_esolver/esolver_lr_lcao_bse.cpp b/source/source_esolver/esolver_lr_lcao_bse.cpp
index be778ecdf20..e0720f003c4 100644
--- a/source/source_esolver/esolver_lr_lcao_bse.cpp
+++ b/source/source_esolver/esolver_lr_lcao_bse.cpp
@@ -516,7 +516,7 @@ void ESolver_BSE<T, TR>::after_all_runners(BaseCell& basecell)
                 for (int is = 0; is < this->X.size(); ++is)
                 {
                     LR::LR_Spectrum<T> spectrum(this->nspin, this->nbasis, this->nocc, this->nvirt, *this->pw_rho, *this->psi_ks,
-                                                *this->ucell_, this->kv, this->gd(), this->orb_cutoff_, this->tcb(),
+                                                *this->ucell_, this->kv, this->kv.wk, this->gd(), this->orb_cutoff_, this->tcb(),
                                                 this->paraX_, this->paraC_, this->paraMat_,
                                                 &this->tda_ene[is * this->nstates], this->eig_ks.c,
                                                 this->X[is].template data<T>(), this->nstates, false/*openshell*/,
@@ -547,7 +547,7 @@ void ESolver_BSE<T, TR>::after_all_runners(BaseCell& basecell)
                 for (int is = 0;is < this->full_X.size();++is)
                 {
                     LR::LR_Spectrum<T> spectrum(this->nspin, this->nbasis, this->nocc, this->nvirt, *this->pw_rho, *this->psi_ks,
-                                                *this->ucell_, this->kv, this->gd(), this->orb_cutoff_, this->tcb(),
+                                                *this->ucell_, this->kv, this->kv.wk, this->gd(), this->orb_cutoff_, this->tcb(),
                                                 this->paraX_, this->paraC_, this->paraMat_,
                                                 &this->full_ene[is * this->nstates], this->eig_ks.c,
                                                 this->full_X[is].template data<T>(), this->nstates, false/*openshell*/,
diff --git a/source/source_esolver/esolver_lr_lcao_tddft.cpp b/source/source_esolver/esolver_lr_lcao_tddft.cpp
index 94427875726..345747232a0 100644
--- a/source/source_esolver/esolver_lr_lcao_tddft.cpp
+++ b/source/source_esolver/esolver_lr_lcao_tddft.cpp
@@ -6,6 +6,7 @@
 #include "source_lcao/module_lr/hamilt_ulr.h"
 #include "source_lcao/module_lr/hsolver_lrtd.h"
 #include "source_lcao/module_lr/lr_spectrum.h"
+#include "source_lcao/module_lr/utils/spectrum_weights.h"
 #include "source_lcao/module_lr/lr_density.h"
 #include "source_lcao/module_lr/ri_benchmark/ri_benchmark.h"
 #include "source_lcao/module_lr/operator_casida/operator_lr_diag.h"
@@ -288,6 +289,7 @@ void ModuleESolver::ESolver_LR<T, TR>::after_all_runners(BaseCell& basecell)
 
     const std::string gauge = LR_Util::tolower(this->inp_->abs_gauge);
     const bool calculate_spectrum = gauge != "none";
+    const std::vector<double> spatial_weights = LR_Util::spatial_k_weights(this->kv.wk, nspin);
     //cal spectrum
     if (LR_Util::tolower(this->inp_->abs_gauge) == "velocity" )
     {
@@ -321,7 +323,7 @@ void ModuleESolver::ESolver_LR<T, TR>::after_all_runners(BaseCell& basecell)
         else
         {
             LR_Spectrum<T> spectrum(nspin, this->nbasis, this->nocc, this->nvirt, *this->pw_rho, *this->psi_ks,
-                *this->ucell_, this->kv, this->gd(), this->orb_cutoff_, this->tcb(),
+                *this->ucell_, this->kv, spatial_weights, this->gd(), this->orb_cutoff_, this->tcb(),
                 this->paraX_, this->paraC_, this->paraMat_,
                 &this->pelec->ekb.c[is * nstates], this->eig_ks.c, this->X[is].template data<T>(), nstates, openshell,
                 LR_Util::tolower(this->inp_->abs_gauge), this->my_rank_, this->out_dir);
diff --git a/source/source_esolver/esolver_lr_orb.cpp b/source/source_esolver/esolver_lr_orb.cpp
index ff0492ddf61..274cd4ca423 100644
--- a/source/source_esolver/esolver_lr_orb.cpp
+++ b/source/source_esolver/esolver_lr_orb.cpp
@@ -4,6 +4,8 @@
 #include "source_lcao/module_lr/ri_benchmark/ri_benchmark.h"
 #include "source_io/module_wf/read_wfc_nao.h"
 #include "source_io/module_output/cube_io.h"
+#include "source_estate/module_charge/chg_rhog_io.h"
+#include "source_base/module_parallel/para_bridge.h"
 #include "source_base/module_external/scalapack_connector.h"
 #include <algorithm>
 
@@ -293,6 +295,22 @@ void ModuleESolver::ESolver_LR<T, TR>::read_ks_chg(Charge& chg_gs)
     chg_gs.set_rhopw(this->pw_rho);
     const bool kin_den = XC_Functional::get_ked_flag() || (this->inp_->out_elf[0] > 0); // mohan add 20251202
     chg_gs.allocate(this->nspin, kin_den, XC_Functional::get_ked_flag(), this->inp_->test_charge);
+    // Keep the existing cube input authoritative. Choose one density source for
+    // all spin channels and broadcast the decision before collective readers.
+    const Parallel::ParaWorld world = Parallel::make_pw_world();
+    const bool have_cubes = module_charge::has_complete_charge_cubes(this->in_dir, this->nspin, world);
+    if (!have_cubes)
+    {
+        const std::string density_file = this->in_dir + this->inp_->suffix + "-CHARGE-DENSITY.restart";
+        const bool loaded = module_charge::read_rho_from_binary(density_file, *this->pw_rho, this->nspin,
+            chg_gs.rhog, chg_gs.rho, world, this->ofs_running_);
+        if (!loaded)
+        {
+            const std::string message = "No complete charge cube set or readable density restart in " + this->in_dir;
+            ModuleBase::WARNING_QUIT("ESolver_LR::read_ks_chg", message);
+        }
+        return;
+    }
     this->ofs_running_ << " try to read charge from file : ";
     for (int is = 0; is < this->nspin; ++is)
     {
diff --git a/source/source_estate/module_charge/chg_rhog_io.cpp b/source/source_estate/module_charge/chg_rhog_io.cpp
index c199db7313a..0e44b165b62 100644
--- a/source/source_estate/module_charge/chg_rhog_io.cpp
+++ b/source/source_estate/module_charge/chg_rhog_io.cpp
@@ -6,6 +6,7 @@
 
 #include <algorithm>
 #include <numeric>
+#include <fstream>
 #include <unistd.h>
 
 namespace
@@ -63,12 +64,15 @@ void module_charge::write_density_restart(Charge& charge,
     }
 }
 
-bool module_charge::read_rhog(const std::string& filename,
-                              const ModulePW::PW_Basis* pw_rhod,
-                              const int nspin,
-                              std::complex<double>** rhog,
-                              const Parallel::ParaWorld& pw_world,
-                              std::ostream* os_warning)
+namespace
+{
+bool read_rhog_impl(const std::string& filename,
+                    const ModulePW::PW_Basis* pw_rhod,
+                    const int nspin,
+                    std::complex<double>** rhog,
+                    const Parallel::ParaWorld& pw_world,
+                    std::ostream* os_warning,
+                    const bool require_same_spin)
 {
     if (pw_rhod == nullptr)
     {
@@ -139,6 +143,12 @@ bool module_charge::read_rhog(const std::string& filename,
         {
             warn(os_warning, pw_world, "module_charge::read_rhog", "some planewaves in file are missing");
         }
+        if ((nspin_in != 1 && nspin_in != 2 && nspin_in != 4) || nspin_in > nspin
+            || (require_same_spin && nspin_in != nspin))
+        {
+            warn(os_warning, pw_world, "module_charge::read_rhog", "unsupported spin channel count in restart");
+            error = true;
+        }
         if (nspin_in < nspin)
         {
             warn(os_warning, pw_world, "module_charge::read_rhog", "some spin channels in file are missing");
@@ -149,7 +159,10 @@ bool module_charge::read_rhog(const std::string& filename,
 
     if (error)
     {
-        warn(os_warning, pw_world, "module_charge::read_rhog", "gamma_only read from file is inconsistent with INPUT");
+        if (gamma_only_in != pw_rhod->gamma_only)
+        {
+            warn(os_warning, pw_world, "module_charge::read_rhog", "gamma_only read from file is inconsistent with INPUT");
+        }
         return false;
     }
 
@@ -230,7 +243,8 @@ bool module_charge::read_rhog(const std::string& filename,
             {
                 int fftixyz = iz + nz * fftixy;
                 int ig = fftixyz2ig[fftixyz];
-                rhog[is][ig] = rhog_in[i];
+                // An owned FFT stick can contain modes outside the current PW cutoff.
+                if (ig >= 0 && ig < pw_rhod->npw) { rhog[is][ig] = rhog_in[i]; }
             }
         }
 
@@ -252,6 +266,59 @@ bool module_charge::read_rhog(const std::string& filename,
     return true;
 }
 
+} // namespace
+
+bool module_charge::read_rhog(const std::string& filename,
+                              const ModulePW::PW_Basis* pw_rhod,
+                              const int nspin,
+                              std::complex<double>** rhog,
+                              const Parallel::ParaWorld& pw_world,
+                              std::ostream* os_warning)
+{
+    return read_rhog_impl(filename, pw_rhod, nspin, rhog, pw_world, os_warning, false);
+}
+
+bool module_charge::has_complete_charge_cubes(const std::string& directory, const int nspin,
+                                             const Parallel::ParaWorld& world)
+{
+    bool complete = nspin == 1 || nspin == 2 || nspin == 4;
+    if (world.rank() == 0)
+    {
+        for (int is = 0; is < nspin && complete; ++is)
+        {
+            const std::string name = nspin == 1 ? "chg.cube" : "chgs" + std::to_string(is + 1) + ".cube";
+            const std::string filename = directory + name;
+            complete = std::ifstream(filename).good();
+        }
+    }
+    Parallel::bcast_bool(complete, world);
+    return complete;
+}
+
+bool module_charge::read_rho_from_binary(const std::string& filename,
+                                      const ModulePW::PW_Basis& rhopw,
+                                      const int nspin,
+                                      std::complex<double>** rhog,
+                                      double** rho,
+                                      const Parallel::ParaWorld& pw_world,
+                                      std::ostream& log)
+{
+    if (rho == nullptr) { return false; }
+    const bool loaded = read_rhog_impl(filename, &rhopw, nspin, rhog, pw_world, &log, true);
+    if (!loaded) { return false; }
+    for (int is = 0; is < nspin; ++is)
+    {
+        rhopw.recip2real(rhog[is], rho[is]);
+    }
+    if (pw_world.rank() == 0)
+    {
+        const std::string message = " Read electron density from binary file and perform inverse FFT: " + filename;
+        log << message << std::endl;
+        std::cout << message << std::endl;
+    }
+    return true;
+}
+
 bool module_charge::write_rhog(const std::string& fchg,
                                const bool gamma_only,
                                const ModulePW::PW_Basis* pw_rho,
diff --git a/source/source_estate/module_charge/chg_rhog_io.h b/source/source_estate/module_charge/chg_rhog_io.h
index e6c2220eb6b..cf0ccdd2b67 100644
--- a/source/source_estate/module_charge/chg_rhog_io.h
+++ b/source/source_estate/module_charge/chg_rhog_io.h
@@ -73,6 +73,23 @@ bool read_rhog(const std::string& filename,
                const Parallel::ParaWorld& pw_world,
                std::ostream* os_warning);
 
+/// Collectively select the complete cube set; partial sets fall back to restart.
+bool has_complete_charge_cubes(const std::string& directory, const int nspin,
+                               const Parallel::ParaWorld& world);
+
+/// Read binary reciprocal-space density and perform inverse FFT for every spin.
+/// Requires the stored spin count to match nspin exactly (no implicit conversion).
+/// Reads density only; does not restore calculation or solver state.
+/// Returns false for recoverable open/configuration failures, leaving rho unchanged.
+/// Malformed or truncated binary input can terminate via Binstream/ WARNING_QUIT.
+bool read_rho_from_binary(const std::string& filename,
+                      const ModulePW::PW_Basis& rhopw,
+                      const int nspin,
+                      std::complex<double>** rhog,
+                      double** rho,
+                      const Parallel::ParaWorld& pw_world,
+                      std::ostream& log);
+
 bool write_rhog(const std::string& fchg,
                 const bool gamma_only,
                 const ModulePW::PW_Basis* pw_rho,
diff --git a/source/source_estate/module_charge/unittests/test_chg_rhog_io.cpp b/source/source_estate/module_charge/unittests/test_chg_rhog_io.cpp
index 3947d1670c2..33bf10c3e8d 100644
--- a/source/source_estate/module_charge/unittests/test_chg_rhog_io.cpp
+++ b/source/source_estate/module_charge/unittests/test_chg_rhog_io.cpp
@@ -4,11 +4,13 @@
 #include "source_base/module_parallel/para_world.h"
 #include "source_base/module_parallel/para_tag.h"
 #include "source_base/module_parallel/para_bridge.h"
+#include "source_base/module_parallel/para_mpi_func.h"
 #ifdef __MPI
 #include "source_basis/module_pw/test/test_tool.h"
 #include "mpi.h"
 #endif
 #include <array>
+#include <cmath>
 #include <complex>
 #include <cstring>
 #include <fstream>
@@ -464,6 +466,212 @@ TEST_F(ReadRhogTest, ReadRhogNspin2To4SpecialPath)
     std::remove(tmpfile.c_str());
 }
 
+TEST_F(ReadRhogTest, RestartDensityTransformsEverySpin)
+{
+    setup_pw_basis();
+    const std::string filename = testing::TempDir() + "lr_density_fft_test.restart";
+    const std::vector<std::array<int, 3>> miller = {{{0, 0, 0}}};
+    const std::vector<std::vector<std::complex<double>>> values = {{{0.25, 0.0}}, {{0.75, 0.0}}};
+    write_binary_by_hand(filename, 0, 2, miller, values);
+    Parallel::barrier(pw_world);
+    fill_rhog(2);
+    std::vector<double> up(rhopw.nrxx, -1.0);
+    std::vector<double> down(rhopw.nrxx, -1.0);
+    double* density[] = {up.data(), down.data()};
+    std::ostringstream log;
+    const bool loaded = module_charge::read_rho_from_binary(filename, rhopw, 2, rhog.data(), density, pw_world, log);
+    EXPECT_TRUE(loaded);
+    for (int ir = 0; ir < rhopw.nrxx; ++ir)
+    {
+        EXPECT_NEAR(up[ir], 0.25, 1e-12);
+        EXPECT_NEAR(down[ir], 0.75, 1e-12);
+    }
+    if (pw_world.rank() == 0)
+    {
+        EXPECT_NE(log.str().find(filename), std::string::npos);
+        std::remove(filename.c_str());
+    }
+}
+
+TEST_F(ReadRhogTest, MissingRestartPreservesRealDensity)
+{
+    setup_pw_basis();
+    fill_rhog(1);
+    const std::string filename = testing::TempDir() + "lr_density_absent_test.restart";
+    if (pw_world.rank() == 0) { std::remove(filename.c_str()); }
+    Parallel::barrier(pw_world);
+    std::vector<double> density(rhopw.nrxx, 0.125);
+    double* channels[] = {density.data()};
+    std::ostringstream log;
+    const bool loaded = module_charge::read_rho_from_binary(filename, rhopw, 1, rhog.data(), channels, pw_world, log);
+    EXPECT_FALSE(loaded);
+    for (double value : density) { EXPECT_DOUBLE_EQ(value, 0.125); }
+}
+
+
+TEST_F(ReadRhogTest, RestartRejectsExcessSpinChannelsBeforeWriting)
+{
+    setup_pw_basis();
+    const std::string filename = testing::TempDir() + "lr_density_spin_mismatch.restart";
+    const std::vector<std::array<int, 3>> miller = {{{0, 0, 0}}};
+    const std::vector<std::vector<std::complex<double>>> values = {{{0.25, 0.0}}, {{0.75, 0.0}}};
+    write_binary_by_hand(filename, 0, 2, miller, values);
+    Parallel::barrier(pw_world);
+    fill_rhog(1); // Exactly one pointer, as allocated by an nspin=1 calculation.
+    const auto original = rhog_data;
+    std::vector<double> density(rhopw.nrxx, 0.125);
+    double* channels[] = {density.data()};
+    std::ostringstream log;
+    EXPECT_FALSE(module_charge::read_rho_from_binary(filename, rhopw, 1, rhog.data(), channels, pw_world, log));
+    EXPECT_EQ(rhog_data, original);
+    for (double value : density) { EXPECT_DOUBLE_EQ(value, 0.125); }
+    Parallel::barrier(pw_world);
+    if (pw_world.rank() == 0) { std::remove(filename.c_str()); }
+}
+
+TEST_F(ReadRhogTest, RestartRejectsMissingSpinChannelsBeforeWriting)
+{
+    setup_pw_basis();
+    const std::string filename = testing::TempDir() + "lr_density_missing_spin.restart";
+    const std::vector<std::array<int, 3>> miller = {{{0, 0, 0}}};
+    const std::vector<std::vector<std::complex<double>>> values = {{{0.25, 0.0}}};
+    write_binary_by_hand(filename, 0, 1, miller, values);
+    Parallel::barrier(pw_world);
+    fill_rhog(2);
+    const auto original = rhog_data;
+    std::vector<double> up(rhopw.nrxx, 0.125);
+    std::vector<double> down(rhopw.nrxx, 0.375);
+    double* channels[] = {up.data(), down.data()};
+    std::ostringstream log;
+    EXPECT_FALSE(module_charge::read_rho_from_binary(filename, rhopw, 2, rhog.data(), channels, pw_world, log));
+    EXPECT_EQ(rhog_data, original);
+    for (double value : up) { EXPECT_DOUBLE_EQ(value, 0.125); }
+    for (double value : down) { EXPECT_DOUBLE_EQ(value, 0.375); }
+    Parallel::barrier(pw_world);
+    if (pw_world.rank() == 0) { std::remove(filename.c_str()); }
+}
+
+TEST_F(ReadRhogTest, RestartSkipsModesOutsideSmallerCutoffOnOwnedStick)
+{
+#ifdef __MPI
+    rhopw.initmpi(pw_world.size(), pw_world.rank(), pw_world.comm());
+#endif
+    rhopw.initgrids(6.5, latvec(), 8, 8, 8);
+    rhopw.initparameters(false, 8.0); // Several sticks for MPI; (0,0,4) lies outside cutoff.
+    rhopw.setuptransform();
+    rhopw.collect_local_pw();
+    ASSERT_GT(rhopw.npwtot, 1);
+    const std::string filename = testing::TempDir() + "lr_density_smaller_cutoff.restart";
+    const std::vector<std::array<int, 3>> miller = {{{0, 0, 0}}, {{0, 0, 1}}, {{0, 0, -1}}, {{0, 0, 4}}};
+    const std::vector<std::vector<std::complex<double>>> values = {{{0.25, 0.0}, {0.1, 0.05}, {0.1, -0.05}, {9.0, 2.0}}};
+    write_binary_by_hand(filename, 0, 1, miller, values);
+    Parallel::barrier(pw_world);
+    // Guard both ends: the former ig=-1 write is visible even without ASan.
+    const std::complex<double> sentinel(77.0, 88.0);
+    std::vector<std::complex<double>> guarded(rhopw.npw + 2, sentinel);
+    std::complex<double>* reciprocal[] = {guarded.data() + 1};
+    std::vector<double> density(rhopw.nrxx, -1.0);
+    double* channels[] = {density.data()};
+    std::ostringstream log;
+    EXPECT_TRUE(module_charge::read_rho_from_binary(filename, rhopw, 1, reciprocal, channels, pw_world, log));
+    EXPECT_EQ(guarded.front(), sentinel);
+    EXPECT_EQ(guarded.back(), sentinel);
+    const double two_pi = 2.0 * std::acos(-1.0);
+    for (int ir = 0; ir < rhopw.nrxx; ++ir)
+    {
+        const int iz = ir % rhopw.nplane + rhopw.startz_current;
+        const double phase = two_pi * iz / rhopw.nz;
+        const double expected = 0.25 + 0.2 * std::cos(phase) - 0.1 * std::sin(phase);
+        EXPECT_NEAR(density[ir], expected, 1e-12);
+    }
+    Parallel::barrier(pw_world);
+    if (pw_world.rank() == 0) { std::remove(filename.c_str()); }
+}
+
+
+TEST_F(ReadRhogTest, CompleteCubeSetTakesPrecedenceOverRestart)
+{
+    const std::string directory = testing::TempDir() + "lr_cube_priority_";
+    const std::string restart = directory + "ABACUS-CHARGE-DENSITY.restart";
+    const std::string cube = directory + "chg.cube";
+    const std::string up_cube = directory + "chgs1.cube";
+    const std::string down_cube = directory + "chgs2.cube";
+    const std::vector<std::array<int, 3>> miller = {{{0, 0, 0}}};
+    const std::vector<std::vector<std::complex<double>>> values = {{{0.25, 0.0}}, {{0.75, 0.0}}};
+    write_binary_by_hand(restart, 0, 2, miller, values);
+    if (pw_world.rank() == 0)
+    {
+        std::ofstream(cube) << "cube fixture";
+        std::ofstream(up_cube) << "cube fixture";
+        std::ofstream(down_cube) << "cube fixture";
+    }
+    Parallel::barrier(pw_world);
+    EXPECT_TRUE(module_charge::has_complete_charge_cubes(directory, 1, pw_world));
+    EXPECT_TRUE(module_charge::has_complete_charge_cubes(directory, 2, pw_world));
+    Parallel::barrier(pw_world);
+    if (pw_world.rank() == 0)
+    {
+        std::remove(cube.c_str());
+        std::remove(up_cube.c_str());
+        std::remove(down_cube.c_str());
+        std::remove(restart.c_str());
+    }
+}
+
+TEST_F(ReadRhogTest, PartialCubeSetUsesRestartForEverySpin)
+{
+    setup_pw_basis();
+    const std::string directory = testing::TempDir() + "lr_cube_partial_";
+    const std::string restart = directory + "ABACUS-CHARGE-DENSITY.restart";
+    const std::string up_cube = directory + "chgs1.cube";
+    const std::string down_cube = directory + "chgs2.cube";
+    const std::vector<std::array<int, 3>> miller = {{{0, 0, 0}}};
+    const std::vector<std::vector<std::complex<double>>> values = {{{0.25, 0.0}}, {{0.75, 0.0}}};
+    write_binary_by_hand(restart, 0, 2, miller, values);
+    if (pw_world.rank() == 0)
+    {
+        std::ofstream(up_cube) << "cube fixture";
+        std::remove(down_cube.c_str());
+    }
+    Parallel::barrier(pw_world);
+    EXPECT_FALSE(module_charge::has_complete_charge_cubes(directory, 2, pw_world));
+    fill_rhog(2);
+    std::vector<double> up(rhopw.nrxx, -1.0);
+    std::vector<double> down(rhopw.nrxx, -1.0);
+    double* density[] = {up.data(), down.data()};
+    std::ostringstream log;
+    EXPECT_TRUE(module_charge::read_rho_from_binary(restart, rhopw, 2, rhog.data(), density, pw_world, log));
+    for (double value : up) { EXPECT_NEAR(value, 0.25, 1e-12); }
+    for (double value : down) { EXPECT_NEAR(value, 0.75, 1e-12); }
+    Parallel::barrier(pw_world);
+    if (pw_world.rank() == 0)
+    {
+        std::remove(up_cube.c_str());
+        std::remove(restart.c_str());
+    }
+}
+
+TEST_F(ReadRhogTest, NeitherDensitySourceIsAvailable)
+{
+    setup_pw_basis();
+    const std::string directory = testing::TempDir() + "lr_cube_neither_";
+    const std::string restart = directory + "ABACUS-CHARGE-DENSITY.restart";
+    const std::string cube = directory + "chg.cube";
+    if (pw_world.rank() == 0)
+    {
+        std::remove(cube.c_str());
+        std::remove(restart.c_str());
+    }
+    Parallel::barrier(pw_world);
+    EXPECT_FALSE(module_charge::has_complete_charge_cubes(directory, 1, pw_world));
+    fill_rhog(1);
+    std::vector<double> rho(rhopw.nrxx, 0.125);
+    double* density[] = {rho.data()};
+    std::ostringstream log;
+    EXPECT_FALSE(module_charge::read_rho_from_binary(restart, rhopw, 1, rhog.data(), density, pw_world, log));
+    for (double value : rho) { EXPECT_DOUBLE_EQ(value, 0.125); }
+}
+
 int main(int argc, char** argv)
 {
 #ifdef __MPI
diff --git a/source/source_io/module_parameter/read_inp_sys.cpp b/source/source_io/module_parameter/read_inp_sys.cpp
index 51c8d95a8b6..f6864a7d35b 100644
--- a/source/source_io/module_parameter/read_inp_sys.cpp
+++ b/source/source_io/module_parameter/read_inp_sys.cpp
@@ -1544,7 +1544,11 @@ updates structures smoothly enough for second-order extrapolation.)";
         item.annotation = "directory of files for reading";
         item.category = "Input files";
         item.type = "String";
-        item.description = "Location of files, such as the electron density (chgs1.cube), required as a starting point.";
+        item.description = "Location of files, such as the electron density (chgs1.cube), required as a starting point. "
+                           "Standalone LR first reads the complete chg.cube/chgs*.cube set. If unavailable, it reads "
+                           "the binary density file {suffix}-CHARGE-DENSITY.restart and performs inverse FFT onto the current grid. "
+                           "The stored spin count must match the requested nspin exactly. "
+                           "This reads ground-state density only; it does not resume a previous calculation.";
         item.default_value = "OUT.$suffix";
         read_sync_string(input.read_file_dir);
         item.reset_value = [](const Input_Item& item, Parameter& para) {
diff --git a/source/source_io/module_restart/restart_exx_csr.cpp b/source/source_io/module_restart/restart_exx_csr.cpp
index 2c7db4a85cd..c565960b834 100644
--- a/source/source_io/module_restart/restart_exx_csr.cpp
+++ b/source/source_io/module_restart/restart_exx_csr.cpp
@@ -21,12 +21,20 @@ namespace ModuleIO
         std::vector<std::map<int, std::map<TAC, RI::Tensor<Tdata>>>>& Hexxs)
     {
         ModuleBase::TITLE("ModuleIO", "read_Hexxs_csr");
+        Hexxs.clear();
         Hexxs.resize(nspin);
         for (int is = 0;is < nspin;++is)
         {
-            ModuleIO::csrFileReader<Tdata> csr(file_name + "_" + std::to_string(is) + ".csr");
+            const std::string path = file_name + "_" + std::to_string(is) + ".csr";
+            ModuleIO::csrFileReader<Tdata> csr(path);
             int nR = csr.getNumberOfR();
-            int nbasis = csr.getMatrixDimension();
+            if (csr.getMatrixDimension() != nbasis)
+            {
+                ModuleBase::WARNING_QUIT("read_Hexxs_csr", "HexxR matrix dimension does not match the LCAO basis");
+            }
+            const std::string confirmation = "Read HexxR from " + path;
+            std::cout << confirmation << std::endl;
+            ModuleBase::TITLE("read_Hexxs_csr", confirmation, false);
             // allocate Hexxs[is]
             for (int iat1 = 0; iat1 < ucell.nat; ++iat1) {
                 for (int iat2 = 0;iat2 < ucell.nat;++iat2) {
diff --git a/source/source_lcao/module_lr/lr_spectrum.h b/source/source_lcao/module_lr/lr_spectrum.h
index 20078a38586..da3df653cce 100644
--- a/source/source_lcao/module_lr/lr_spectrum.h
+++ b/source/source_lcao/module_lr/lr_spectrum.h
@@ -17,14 +17,14 @@ namespace LR
     public:
         LR_Spectrum(const int& nspin_global, const int& naos, const std::vector<int>& nocc, const std::vector<int>& nvirt,
             const ModulePW::PW_Basis& rho_basis, psi::Psi<T>& psi_ks_in,
-            const UnitCell& ucell, const K_Vectors& kv_in, const Grid_Driver& gd, const std::vector<double>& orb_cutoff,
+            const UnitCell& ucell, const K_Vectors& kv_in, const std::vector<double>& spatial_weights_in, const Grid_Driver& gd, const std::vector<double>& orb_cutoff,
             const TwoCenterBundle& two_center_bundle_,
             const std::vector<Parallel_2D>& pX_in, const Parallel_2D& pc_in, const Parallel_Orbitals& pmat_in,
             const double* omega, const double* eig_ks, const T* X, const int& nstate, const bool& openshell,
             const std::string& gauge, const int my_rank, const std::string& out_dir) :
             my_rank(my_rank), nspin_x(openshell ? 2 : 1), naos(naos), nocc(nocc), nvirt(nvirt),
             nk(nspin_global == 2 ? kv_in.get_nks() / 2 : kv_in.get_nks()),
-            rho_basis(rho_basis), ucell(ucell), kv(kv_in), gd_(gd),
+            rho_basis(rho_basis), ucell(ucell), kv(kv_in), spatial_weights(spatial_weights_in), gd_(gd),
             orb_cutoff_(orb_cutoff), two_center_bundle_(two_center_bundle_),
             pX(pX_in), pc(pc_in), pmat(pmat_in),
             omega(omega), eig_ks(eig_ks), X(X), nstate(nstate), gauge(gauge), out_dir(out_dir),
@@ -92,6 +92,8 @@ namespace LR
         T* Y = nullptr; ///< the deexcitation part of amplitudes
         bool is_full = false;
         const K_Vectors& kv;
+        /// Pure spatial weights, with no spin degeneracy; same spin/k ordering as kv.
+        const std::vector<double> spatial_weights;
         std::vector<psi::Psi<T>> psi_ks_vec;
         const std::vector<Parallel_2D>& pX;
         const Parallel_2D& pc;
diff --git a/source/source_lcao/module_lr/lr_spectrum_velocity.cpp b/source/source_lcao/module_lr/lr_spectrum_velocity.cpp
index 96a0a4612e4..902b73e7697 100644
--- a/source/source_lcao/module_lr/lr_spectrum_velocity.cpp
+++ b/source/source_lcao/module_lr/lr_spectrum_velocity.cpp
@@ -1,4 +1,5 @@
 #include "lr_spectrum.h"
+#include "utils/spectrum_weights.h"
 #include "source_base/global_variable.h"
 #include "source_lcao/module_lr/dm_trans/dm_trans.h"
 #include "source_lcao/module_lr/utils/lr_util_hcontainer.h"
@@ -161,7 +162,6 @@ namespace LR
                     const int goffset_v_ds = (is * 3 + id) * nk * nbands * nbands;
                     for (int ik = 0; ik < nk; ++ik)
                     {
-                        const double wk = this->kv.wk[ik] * nk; // k-point weight, normalized as sum = nk
                         const std::size_t loffset_X = loffset_X_bs + ik * pX[is].get_local_size();
                         const int goffset_v = goffset_v_ds + ik * nbands * nbands;
                         for (int io = 0; io < pX[is].get_col_size(); ++io)    // nocc_local
@@ -172,29 +172,10 @@ namespace LR
                                 int iv_g = pX[is].local2global_row(iv);
                                 const std::size_t X_index = loffset_X + io * pX[is].get_row_size() + iv;
                                 const int v_index = goffset_v + (iv_g+nocc[is]) * nbands + io_g;
-                                if (use_ks_gap)
-                                {
-                                    td += this->vmo_ptr[v_index] * X[X_index] / eig_ks_diff[X_index - loffset_X_b] * wk;
-                                    if (this->is_full)
-                                    {
-                                        // THE HERMITIAN CONJUGATE OF VMO HAS BEEN VERIFIED
-                                        // const int v_index2 = goffset_v + io_g * nbands + iv_g + nocc[is];
-                                        // if (std::abs(vmo_ptr[v_index2] - std::conj(this->vmo_ptr[v_index])) > 1e-5){
-                                        //     std::cout<<"io:"<<io_g<<" iv:"<<iv_g<<" v:"<<vmo_ptr[v_index]
-                                        //         <<" v^T:" << vmo_ptr[v_index2]<<std::endl;
-                                        //     }
-                                        // <ik|v|ak>X_{aik}/(Ea-Ei) + <ak|v|ik>Y_{aik}/(Ei-Ea)
-                                        td -= std::conj(this->vmo_ptr[v_index]) * Y[X_index] / eig_ks_diff[X_index - loffset_X_b] * wk;
-                                    }
-                                }
-                                else
-                                {
-                                    td += this->vmo_ptr[v_index] * X[X_index] * wk;
-                                    if (this->is_full)
-                                    {
-                                        td += std::conj(this->vmo_ptr[v_index]) * Y[X_index] * wk;
-                                    }
-                                }
+                                const T deexcitation = this->is_full ? Y[X_index] : T(0);
+                                const double gap = use_ks_gap ? eig_ks_diff[X_index - loffset_X_b] : 1.0;
+                                td += LR_Util::weighted_velocity_transition(this->vmo_ptr[v_index], X[X_index],
+                                    deexcitation, gap, use_ks_gap, this->spatial_weights, is, ik, nk);
                             }
                         }
                     }
diff --git a/source/source_lcao/module_lr/test/CMakeLists.txt b/source/source_lcao/module_lr/test/CMakeLists.txt
index 3c019fe7df2..c6180ff87ed 100644
--- a/source/source_lcao/module_lr/test/CMakeLists.txt
+++ b/source/source_lcao/module_lr/test/CMakeLists.txt
@@ -161,3 +161,8 @@ if(ENABLE_MPI)
     PROPERTIES PROCESSORS 4 TIMEOUT 60 ENVIRONMENT
                "OMP_NUM_THREADS=1;MKL_NUM_THREADS=1;OPENBLAS_NUM_THREADS=1")
 endif()
+
+AddTest(
+  TARGET MODULE_LR_spectrum_weights
+  SOURCES test_spectrum_weights.cpp
+)
diff --git a/source/source_lcao/module_lr/test/test_spectrum_weights.cpp b/source/source_lcao/module_lr/test/test_spectrum_weights.cpp
new file mode 100644
index 00000000000..a15453934f2
--- /dev/null
+++ b/source/source_lcao/module_lr/test/test_spectrum_weights.cpp
@@ -0,0 +1,105 @@
+#include "../utils/spectrum_weights.h"
+#include <gtest/gtest.h>
+#include <cmath>
+#include <complex>
+
+TEST(SpectrumWeights, RemovesOnlyUnpolarizedSpinDegeneracy)
+{
+    EXPECT_DOUBLE_EQ(LR_Util::spatial_k_weight(0.002, 1), 0.001);
+    EXPECT_DOUBLE_EQ(LR_Util::spatial_k_weight(0.001, 2), 0.001);
+    EXPECT_DOUBLE_EQ(LR_Util::spatial_k_weight(0.001, 4), 0.001);
+}
+
+TEST(SpectrumWeights, PreservesNonuniformSpatialWeights)
+{
+    const double ks_weights[] = {0.2, 0.6, 1.2};
+    const double expected[] = {0.1, 0.3, 0.6};
+    double sum = 0.0;
+    for (int k = 0; k < 3; ++k)
+    {
+        const double w = LR_Util::spatial_k_weight(ks_weights[k], 1);
+        EXPECT_DOUBLE_EQ(w, expected[k]);
+        sum += w;
+    }
+    EXPECT_DOUBLE_EQ(sum, 1.0);
+}
+
+TEST(SpectrumWeights, SingletIntensityMatchesExplicitSpinChannels)
+{
+    // A unit-normalized singlet has X_up=X_down=X/sqrt(2).
+    // Compare that explicit contraction with the spatial-orbital sqrt(2) path.
+    const int nk = 2;
+    const std::complex<double> dipoles[] = {{1.0, 0.5}, {-0.2, 0.7}};
+    const double amplitudes[] = {0.6, 0.8};
+    std::complex<double> singlet = 0.0;
+    std::complex<double> explicit_spins = 0.0;
+    for (int k = 0; k < nk; ++k)
+    {
+        const double unpolarized_w = LR_Util::spatial_k_weight(1.0, 1);
+        const double polarized_w = LR_Util::spatial_k_weight(0.5, 2);
+        singlet += std::sqrt(2.0) * dipoles[k] * amplitudes[k] * unpolarized_w * double(nk);
+        for (int spin = 0; spin < 2; ++spin)
+        {
+            explicit_spins += dipoles[k] * amplitudes[k] / std::sqrt(2.0) * polarized_w * double(nk);
+        }
+    }
+    EXPECT_NEAR(std::abs(singlet - explicit_spins), 0.0, 1e-14);
+    EXPECT_NEAR(std::norm(singlet), 1.8664, 1e-14);
+}
+
+TEST(SpectrumWeights, LRConversionPreservesSharedKSWeights)
+{
+    const std::vector<double> ks_weights = {0.2, 0.6, 1.2};
+    const std::vector<double> spatial = LR_Util::spatial_k_weights(ks_weights, 1);
+    EXPECT_EQ(spatial, (std::vector<double>{0.1, 0.3, 0.6}));
+    EXPECT_EQ(ks_weights, (std::vector<double>{0.2, 0.6, 1.2}));
+    EXPECT_EQ(LR_Util::spatial_k_weights(spatial, 2), spatial);
+}
+
+TEST(SpectrumWeights, LRAndBSEIndependentTransitionHaveSameSingletStrength)
+{
+    // BSE RI_kRlist supplies 1/Nk directly; KS supplies 2/Nk for nspin=1.
+    // An IPA eigenvector has one unit entry, with no compensating spin scale.
+    for (const int nk : {1, 2, 8})
+    {
+        const std::vector<double> ks_weights(nk, 2.0 / nk);
+        const std::vector<double> lr_weights = LR_Util::spatial_k_weights(ks_weights, 1);
+        const std::vector<double> bse_weights(nk, 1.0 / nk);
+        const std::complex<double> dipole(1.0, 0.5);
+        const std::complex<double> lr_amplitude = std::sqrt(2.0) * LR_Util::weighted_velocity_transition(
+            dipole, 1.0, 0.0, 1.0, true, lr_weights, 0, 0, nk);
+        const std::complex<double> bse_amplitude = std::sqrt(2.0) * LR_Util::weighted_velocity_transition(
+            dipole, 1.0, 0.0, 1.0, true, bse_weights, 0, 0, nk);
+        EXPECT_NEAR(std::norm(lr_amplitude), 2.5, 1e-14);
+        EXPECT_NEAR(std::abs(lr_amplitude - bse_amplitude), 0.0, 1e-14);
+    }
+}
+
+TEST(SpectrumWeights, ProductionContractionSelectsEachSpinBlock)
+{
+    const std::vector<double> ks_weights = {0.1, 0.9, 0.3, 0.7};
+    const std::vector<double> weights = LR_Util::spatial_k_weights(ks_weights, 2);
+    const std::complex<double> velocity(2.0, 4.0);
+    const double expected[] = {0.2, 1.8, 0.6, 1.4};
+    for (int spin = 0; spin < 2; ++spin)
+    {
+        for (int ik = 0; ik < 2; ++ik)
+        {
+            const std::complex<double> actual = LR_Util::weighted_velocity_transition(
+                velocity, 1.0, 0.0, 2.0, true, weights, spin, ik, 2);
+            EXPECT_NEAR(std::abs(actual - std::complex<double>(1.0, 2.0) * expected[spin * 2 + ik]), 0.0, 1e-14);
+        }
+    }
+}
+
+TEST(SpectrumWeights, ProductionFullBSEContractionKeepsDeexcitationSign)
+{
+    const std::vector<double> bse_weights = {1.0};
+    const std::complex<double> velocity(2.0, 4.0);
+    const std::complex<double> ks_gap = LR_Util::weighted_velocity_transition(
+        velocity, 3.0, 1.0, 2.0, true, bse_weights, 0, 0, 1);
+    const std::complex<double> lr_energy = LR_Util::weighted_velocity_transition(
+        velocity, 3.0, 1.0, 2.0, false, bse_weights, 0, 0, 1);
+    EXPECT_EQ(ks_gap, std::complex<double>(2.0, 8.0));
+    EXPECT_EQ(lr_energy, std::complex<double>(8.0, 8.0));
+}
diff --git a/source/source_lcao/module_lr/utils/spectrum_weights.h b/source/source_lcao/module_lr/utils/spectrum_weights.h
new file mode 100644
index 00000000000..0464f53d17f
--- /dev/null
+++ b/source/source_lcao/module_lr/utils/spectrum_weights.h
@@ -0,0 +1,48 @@
+#ifndef ABACUS_LR_SPECTRUM_WEIGHTS_H
+#define ABACUS_LR_SPECTRUM_WEIGHTS_H
+
+#include <vector>
+#include <complex>
+
+namespace LR_Util
+{
+/// KS weights include a factor two only for nspin=1. Optical singlet
+/// amplitudes account for that degeneracy separately through sqrt(2).
+inline double spatial_k_weight(const double ks_weight, const int nspin_global)
+{
+    return nspin_global == 1 ? 0.5 * ks_weight : ks_weight;
+}
+
+/// Convert KS weights at the LR/spectrum boundary without changing the shared k list.
+/// BSE already owns spatial weights and must pass them directly to LR_Spectrum.
+inline std::vector<double> spatial_k_weights(const std::vector<double>& ks_weights, const int nspin_global)
+{
+    std::vector<double> weights(ks_weights);
+    for (double& weight : weights)
+    {
+        weight = spatial_k_weight(weight, nspin_global);
+    }
+    return weights;
+}
+
+/// One production velocity-gauge transition, before the common i/omega and
+/// singlet factors. Spatial weights are indexed by response spin and k point.
+template <typename T>
+inline std::complex<double> weighted_velocity_transition(
+    const std::complex<double>& velocity, const T& x, const T& y,
+    const double gap, const bool use_ks_gap, const std::vector<double>& weights,
+    const int spin, const int ik, const int nk)
+{
+    const double weight = weights[spin * nk + ik] * nk;
+    const std::complex<double> excitation = velocity * x;
+    const std::complex<double> deexcitation = std::conj(velocity) * y;
+    if (use_ks_gap)
+    {
+        return (excitation - deexcitation) * (weight / gap);
+    }
+    return (excitation + deexcitation) * weight;
+}
+
+}
+
+#endif