Write DFT+U occupation eigenvalues to occ_mat.txt

source/source_pw/module_pwdft/dftu_base_io.cpp

@@ -293,6 +293,7 @@ void write_latest_occmat(const Plus_U_Base& dftu,
     // esolver gates on latest_ready, the LCAO esolver computes it every
     // electronic iteration.
     write_snapshot_section(ofdftu, dftu, ucell, step, true, soc_layout);
+    write_occup_eigenvalues(dftu, ucell, ofdftu, soc_layout);
 }
 
 

source/source_pw/module_pwdft/dftu_base_io.h

@@ -149,6 +149,12 @@ void write_occup_m(const Plus_U_Base& dftu,
                    std::ostream& ofs,
                    OccmatSocLayout soc_layout);
 
+/// Append sorted occupation eigenvalues as restart-safe comment lines.
+void write_occup_eigenvalues(const Plus_U_Base& dftu,
+                             const UnitCell& ucell,
+                             std::ostream& ofs,
+                             OccmatSocLayout soc_layout);
+
 } // namespace DFTU_BASE
 
 #endif

source/source_pw/module_pwdft/dftu_base_write.cpp

@@ -23,10 +23,7 @@ struct SocCompactBlocks
     std::vector<double> down_down;
 };
 
-// Retained for the follow-up occupation-eigensystem output refactor. The
-// compact writer does not currently call these helpers: this PR only removes
-// the legacy text format and separates the writer from the file orchestration.
-namespace retained_eigenvalue_helpers
+namespace
 {
 
 inline void jacobi_rotate(std::vector<std::vector<double>>& matrix,
@@ -67,6 +64,11 @@ inline std::vector<double> calculate_eigenvalues(std::vector<std::vector<double>
                                                  const int size)
 {
     std::vector<double> eigenvalues(size);
+    if (size == 1)
+    {
+        eigenvalues[0] = matrix[0][0];
+        return eigenvalues;
+    }
     while (true)
     {
         int p = 0;
@@ -93,6 +95,7 @@ inline std::vector<double> calculate_eigenvalues(std::vector<std::vector<double>
         }
         jacobi_rotate(matrix, p, q, size);
     }
+    std::sort(eigenvalues.begin(), eigenvalues.end());
     return eigenvalues;
 }
 
@@ -159,7 +162,7 @@ inline std::vector<double> calculate_spinor_eigenvalues(const SocCompactBlocks&
     return eigenvalues;
 }
 
-} // namespace retained_eigenvalue_helpers
+} // namespace
 
 static SocCompactBlocks extract_soc_compact_blocks(const Plus_U_Base& dftu,
                                                    const int iat,
@@ -326,4 +329,62 @@ void write_occup_m(const Plus_U_Base& dftu,
     }
 }
 
+void write_occup_eigenvalues(const Plus_U_Base& dftu,
+                             const UnitCell& ucell,
+                             std::ostream& os,
+                             const OccmatSocLayout soc_layout)
+{
+    const int nspin = dftu.occmat().nspin();
+    for (int type = 0; type < ucell.ntype; ++type)
+    {
+        if (!dftu.has_l_channel(type))
+        {
+            continue;
+        }
+        const int l = dftu.get_l_channel(type);
+        const int size = 2 * l + 1;
+        for (int atom = 0; atom < ucell.atoms[type].na; ++atom)
+        {
+            const int iat = ucell.itia2iat(type, atom);
+            const int channels = nspin == 4 ? 1 : 2;
+            for (int spin = 0; spin < channels; ++spin)
+            {
+                std::vector<double> values;
+                if (nspin == 4)
+                {
+                    const SocCompactBlocks blocks = extract_soc_compact_blocks(dftu, iat, l, soc_layout);
+                    values = calculate_spinor_eigenvalues(blocks, size);
+                }
+                else
+                {
+                    std::vector<std::vector<double>> matrix(size, std::vector<double>(size));
+                    for (int row = 0; row < size; ++row)
+                    {
+                        for (int column = 0; column < size; ++column)
+                        {
+                            matrix[row][column] = dftu.occmat().get(iat, l, spin, row, column);
+                        }
+                    }
+                    values = calculate_eigenvalues(matrix, size);
+                }
+                os << "# eigenvalues atom " << iat + 1 << " L " << l;
+                if (nspin == 4)
+                {
+                    os << " spinor:";
+                }
+                else
+                {
+                    os << " spin " << spin + 1 << ":";
+                }
+                os << std::fixed << std::setprecision(8);
+                for (const double value : values)
+                {
+                    os << " " << value;
+                }
+                os << "\n";
+            }
+        }
+    }
+}
+
 } // namespace DFTU_BASE

source/source_pw/module_pwdft/unittests/test_dftu_base.cpp

@@ -19,6 +19,7 @@
 #endif
 
 #include <cstdio>
+#include <cmath>
 #include <fstream>
 #include <sstream>
 #include <string>
@@ -530,6 +531,7 @@ TEST_F(DFTUBaseTest, AppendSnapshotNAPlaceholderAndReady)
     const std::string content2 = ss2.str();
     EXPECT_NE(content2.find("Fe Atom 1 L 2"), std::string::npos);
     EXPECT_EQ(content2.find("N/A"), std::string::npos);
+    EXPECT_EQ(content2.find("# eigenvalues"), std::string::npos);
 
     // Clean up so subsequent tests checking "file not written" are not
     // fooled by leftovers.  std::remove is Rule 17 compliant.
@@ -913,6 +915,89 @@ TEST_F(DFTUBaseTest, CompactWriterReaderCollinearRoundtrip)
     EXPECT_NEAR(parsed_nelec[1], expected_nelec[1], 1e-7);
 }
 
+TEST_F(DFTUBaseTest, LatestOccMatIncludesSortedCollinearEigenvalues)
+{
+    Plus_U_Base dftu;
+    init_dftu(dftu, false);
+    for (int spin = 0; spin < 2; ++spin)
+    {
+        for (int m = 0; m < 5; ++m)
+        {
+            dftu.occmat().set(0, 2, spin, m, m, spin == 0 ? 0.8 : 0.2);
+        }
+    }
+    dftu.occmat().set(0, 2, 0, 0, 1, 0.1);
+    dftu.occmat().set(0, 2, 0, 1, 0, 0.1);
+
+    const DFTU_BASE::OccmatOutputCfg cfg = {0, 1, 5, true, 1};
+    const DFTU_BASE::OccmatStepInfo step = {0, 2, 1e-6, 1e-7};
+    const std::string filename = testing::TempDir() + "occ_mat.txt";
+    DFTU_BASE::write_latest_occmat(dftu, ucell, testing::TempDir(), step, cfg,
+                                   DFTU_BASE::SOC_LAYOUT_PAULI);
+
+    std::ifstream input(filename.c_str());
+    ASSERT_TRUE(input.is_open());
+    std::string line;
+    int eigenvalue_lines = 0;
+    while (std::getline(input, line))
+    {
+        if (line.find("# eigenvalues atom 1 L 2 spin ") != 0)
+        {
+            continue;
+        }
+        ++eigenvalue_lines;
+        ASSERT_LE(eigenvalue_lines, 2);
+        std::istringstream values(line.substr(line.find(':') + 1));
+        const double expected[2][5] = {{0.7, 0.8, 0.8, 0.8, 0.9},
+                                       {0.2, 0.2, 0.2, 0.2, 0.2}};
+        double value = 0.0;
+        for (int m = 0; m < 5; ++m)
+        {
+            ASSERT_TRUE(values >> value);
+            EXPECT_NEAR(value, expected[eigenvalue_lines - 1][m], 1e-8);
+        }
+        EXPECT_FALSE(values >> value);
+    }
+    EXPECT_EQ(eigenvalue_lines, 2);
+
+    OccupationMatrix restored;
+    const std::vector<int> l_channel = {2};
+    restored.init(ucell, l_channel, 2, 1);
+    DFTU_BASE::read_occup_m(ucell, restored, l_channel, filename, 2, 1,
+                            DFTU_BASE::SOC_LAYOUT_PAULI);
+    EXPECT_NEAR(restored.get(0, 2, 0, 0, 1), 0.1, 1e-8);
+    std::remove(filename.c_str());
+}
+
+TEST_F(DFTUBaseTest, SocEigenvaluesIncludeImaginarySpinCoupling)
+{
+    Plus_U_Base dftu;
+    init_dftu_for_spin(dftu, false, 4, 2, DFTU_BASE::SOC_LAYOUT_PAULI);
+    const int block_size = 25;
+    for (int m = 0; m < 5; ++m)
+    {
+        const int index = m * 5 + m;
+        dftu.occmat().mat(0, 2, 0).c[index] = 1.0;
+        dftu.occmat().mat(0, 2, 0).c[2 * block_size + index] = 0.6;
+        dftu.occmat().mat(0, 2, 0).c[3 * block_size + index] = 0.6;
+    }
+    std::ostringstream output;
+    DFTU_BASE::write_occup_eigenvalues(dftu, ucell, output,
+                                       DFTU_BASE::SOC_LAYOUT_PAULI);
+    const std::string line = output.str();
+    ASSERT_EQ(line.find("# eigenvalues atom 1 L 2 spinor:"), 0);
+    std::istringstream values(line.substr(line.find(':') + 1));
+    const double low = 0.5 - std::sqrt(0.18);
+    const double high = 0.5 + std::sqrt(0.18);
+    double value = 0.0;
+    for (int i = 0; i < 10; ++i)
+    {
+        ASSERT_TRUE(values >> value);
+        EXPECT_NEAR(value, i < 5 ? low : high, 1e-8);
+    }
+    EXPECT_FALSE(values >> value);
+}
+
 TEST_F(DFTUBaseTest, CompactWriterReaderSocPauliRoundtrip)
 {
     Plus_U_Base source;