From d1fba2e0e5e5e560320ca79a11bbb873f0acad4b Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Wed, 4 Nov 2026 20:08:02 +0800 Subject: [PATCH] Write DFT+U occupation eigenvalues to occ_mat.txt --- .../source_pw/module_pwdft/dftu_base_io.cpp | 1 + source/source_pw/module_pwdft/dftu_base_io.h | 6 ++ .../module_pwdft/dftu_base_write.cpp | 71 ++++++++++++++-- .../module_pwdft/unittests/test_dftu_base.cpp | 85 +++++++++++++++++++ 4 files changed, 158 insertions(+), 5 deletions(-) diff --git a/source/source_pw/module_pwdft/dftu_base_io.cpp b/source/source_pw/module_pwdft/dftu_base_io.cpp index fa8023dfaf6..909ee029fe3 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.cpp +++ b/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); } diff --git a/source/source_pw/module_pwdft/dftu_base_io.h b/source/source_pw/module_pwdft/dftu_base_io.h index 4750cf0fbe7..4de4e6d64e8 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.h +++ b/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 diff --git a/source/source_pw/module_pwdft/dftu_base_write.cpp b/source/source_pw/module_pwdft/dftu_base_write.cpp index d10fe5038c5..10d90680253 100644 --- a/source/source_pw/module_pwdft/dftu_base_write.cpp +++ b/source/source_pw/module_pwdft/dftu_base_write.cpp @@ -23,10 +23,7 @@ struct SocCompactBlocks std::vector 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>& matrix, @@ -67,6 +64,11 @@ inline std::vector calculate_eigenvalues(std::vector const int size) { std::vector eigenvalues(size); + if (size == 1) + { + eigenvalues[0] = matrix[0][0]; + return eigenvalues; + } while (true) { int p = 0; @@ -93,6 +95,7 @@ inline std::vector calculate_eigenvalues(std::vector } jacobi_rotate(matrix, p, q, size); } + std::sort(eigenvalues.begin(), eigenvalues.end()); return eigenvalues; } @@ -159,7 +162,7 @@ inline std::vector 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 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> matrix(size, std::vector(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 diff --git a/source/source_pw/module_pwdft/unittests/test_dftu_base.cpp b/source/source_pw/module_pwdft/unittests/test_dftu_base.cpp index 233af229513..1d4a9ccc60e 100644 --- a/source/source_pw/module_pwdft/unittests/test_dftu_base.cpp +++ b/source/source_pw/module_pwdft/unittests/test_dftu_base.cpp @@ -19,6 +19,7 @@ #endif #include +#include #include #include #include @@ -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 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;