From 2e019ee16407396bd476702d19a21ab3577ec185 Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Wed, 21 Oct 2026 20:29:26 +0800 Subject: [PATCH 1/6] Test DFT+U compact occupation writer roundtrips --- .../module_pwdft/test/dftu_base_test.cpp | 121 +++++++++++++++++- 1 file changed, 118 insertions(+), 3 deletions(-) diff --git a/source/source_pw/module_pwdft/test/dftu_base_test.cpp b/source/source_pw/module_pwdft/test/dftu_base_test.cpp index 56f80a8f04..e5197b98eb 100644 --- a/source/source_pw/module_pwdft/test/dftu_base_test.cpp +++ b/source/source_pw/module_pwdft/test/dftu_base_test.cpp @@ -86,12 +86,25 @@ class DFTUBaseTest : public testing::Test /// Call init_base with the given Yukawa switch on a fresh d orbital void init_dftu(Plus_U_Base& dftu, const bool yukawa_potential) + { + init_dftu_for_spin(dftu, + yukawa_potential, + 2, + 1, + DFTU_BASE::SOC_LAYOUT_PAULI); + } + + void init_dftu_for_spin(Plus_U_Base& dftu, + const bool yukawa_potential, + const int nspin, + const int npol, + const DFTU_BASE::OccmatSocLayout soc_layout) { const std::vector l_channel = {2}; const std::vector hubbard_u = {0.0}; dftu.init_base(ucell, - 1, // npol - 2, // nspin + npol, + nspin, l_channel, yukawa_potential, 0.5, // yukawa_lambda @@ -103,7 +116,7 @@ class DFTUBaseTest : public testing::Test 0.0, // uramping 0, // init_occ_mat 0, // mixing_dftu - DFTU_BASE::SOC_LAYOUT_PAULI); + soc_layout); } }; @@ -829,6 +842,108 @@ TEST_F(DFTUBaseTest, InitBaseReadsOccMatSocPauliRoundtrip) } } +TEST_F(DFTUBaseTest, CompactWriterReaderCollinearRoundtrip) +{ + Plus_U_Base source; + init_dftu_for_spin(source, false, 2, 1, DFTU_BASE::SOC_LAYOUT_PAULI); + for (int spin = 0; spin < 2; ++spin) + { + for (int row = 0; row < 5; ++row) + { + for (int column = 0; column < 5; ++column) + { + const double value = 0.1 * (spin + 1) + 0.01 * row + 0.001 * column; + source.occmat().set(0, 2, spin, row, column, value); + } + } + } + + const std::string filename = testing::TempDir() + "compact_writer_collinear.txt"; + { + std::ofstream output(filename.c_str()); + ASSERT_TRUE(output.is_open()); + DFTU_BASE::write_occup_m(source, + ucell, + output, + true, + 2, + 1, + DFTU_BASE::OCMAT_FMT_READABLE, + DFTU_BASE::SOC_LAYOUT_PAULI); + } + + 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); + + for (int spin = 0; spin < 2; ++spin) + { + for (int row = 0; row < 5; ++row) + { + for (int column = 0; column < 5; ++column) + { + EXPECT_NEAR(restored.get(0, 2, spin, row, column), + source.occmat().get(0, 2, spin, row, column), + 1e-8); + } + } + } +} + +TEST_F(DFTUBaseTest, CompactWriterReaderSocPauliRoundtrip) +{ + Plus_U_Base source; + init_dftu_for_spin(source, false, 4, 2, DFTU_BASE::SOC_LAYOUT_PAULI); + const int block_size = 25; + for (int block = 0; block < 4; ++block) + { + for (int index = 0; index < block_size; ++index) + { + source.occmat().mat(0, 2, 0).c[block * block_size + index] + = 0.1 * (block + 1) + 0.001 * index; + } + } + + const std::string filename = testing::TempDir() + "compact_writer_soc_pauli.txt"; + { + std::ofstream output(filename.c_str()); + ASSERT_TRUE(output.is_open()); + DFTU_BASE::write_occup_m(source, + ucell, + output, + true, + 4, + 2, + DFTU_BASE::OCMAT_FMT_READABLE, + DFTU_BASE::SOC_LAYOUT_PAULI); + } + + OccupationMatrix restored; + const std::vector l_channel = {2}; + restored.init(ucell, l_channel, 4, 2); + DFTU_BASE::read_occup_m(ucell, + restored, + l_channel, + filename, + 4, + 2, + DFTU_BASE::SOC_LAYOUT_PAULI); + + for (int index = 0; index < 4 * block_size; ++index) + { + EXPECT_NEAR(restored.mat(0, 2, 0).c[index], + source.occmat().mat(0, 2, 0).c[index], + 1e-8); + } +} + TEST_F(DFTUBaseTest, ReadOccMatCompactHandlesLongComment) { const std::string fn = testing::TempDir() + "compact_long_comment.txt"; From eaf3c24a8f7acac91e786928c68fa3246591bd86 Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Wed, 21 Oct 2026 22:47:06 +0800 Subject: [PATCH 2/6] Remove legacy DFT+U occupation writer --- .../source_pw/module_pwdft/dftu_base_io.cpp | 337 ++---------------- source/source_pw/module_pwdft/dftu_base_io.h | 12 +- source/source_pw/module_pwdft/dftu_types.h | 7 - .../module_pwdft/test/dftu_base_test.cpp | 6 - 4 files changed, 40 insertions(+), 322 deletions(-) diff --git a/source/source_pw/module_pwdft/dftu_base_io.cpp b/source/source_pw/module_pwdft/dftu_base_io.cpp index a8f680b33e..db5ffca1fc 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.cpp +++ b/source/source_pw/module_pwdft/dftu_base_io.cpp @@ -2,7 +2,6 @@ #include "source_cell/unitcell.h" #include "source_pw/module_pwdft/dftu_base.h" -#include "source_base/constants.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" #include "source_base/module_external/lapack_connector.h" @@ -20,11 +19,12 @@ #include #include -// local helpers for eigenvalue calculation -// migrated from dftu_base.cpp, mohan 2025-11-08 namespace { +// 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 file orchestration. inline void JacobiRotate(std::vector>& A, int p, int q, int n) { if (std::abs(A[p][q]) > 1e-10) @@ -147,16 +147,6 @@ inline void extract_soc_pauli_blocks(const Plus_U_Base& dftu, } } -/// @brief Diagonalize the SOC occupation matrix reconstructed in spinor space. -/// -/// The Hermitian spinor matrix is built from the Pauli blocks as -/// n_upup = (b0 + b3)/2, n_dndn = (b0 - b3)/2, -/// n_updown = (b1 + i*b2)/2. -/// Its 2*(2l+1) real eigenvalues are computed with LAPACK zheev. -/// -/// @param blocks four Pauli-component blocks -/// @param m 2l+1 -/// @return eigenvalues sorted ascending by zheev inline std::vector calculate_spinor_eigenvalues( const std::vector>& blocks, const int m) @@ -434,8 +424,7 @@ void write_snapshot_section(std::ostream& os, } // Full occupation matrices together with per-atom magnetism. - write_occup_m(dftu, ucell, os, true, nspin, npol, - OCMAT_FMT_READABLE, soc_layout); + write_occup_m(dftu, ucell, os, nspin, soc_layout); } void append_ion_step_snapshot(const Plus_U_Base& dftu, @@ -581,10 +570,7 @@ void write_latest_occmat(const Plus_U_Base& dftu, void write_occup_m(const Plus_U_Base& dftu, const UnitCell& ucell, std::ostream& ofs, - bool diag, int nspin, - int npol, - OccmatTextFormat fmt, OccmatSocLayout soc_layout) { ModuleBase::TITLE("DFTU_BASE", "write_occup_m"); @@ -600,170 +586,53 @@ void write_occup_m(const Plus_U_Base& dftu, { continue; } - const int NL = ucell.atoms[T].nwl + 1; - const int LC = dftu.get_l_channel(T); + const int l = dftu.get_l_channel(T); for (int I = 0; I < ucell.atoms[T].na; I++) { const int iat = ucell.itia2iat(T, I); + const std::string& elem_label = ucell.atoms[T].label; - for (int l = 0; l < NL; l++) + if (nspin == 1 || nspin == 2) { - if (l != dftu.get_l_channel(T)) - { - continue; - } - - const std::string& elem_label = ucell.atoms[T].label; - const int index_in_type = I + 1; - const ModuleBase::Vector3& tau = ucell.atoms[T].tau[I]; - const bool readable = (fmt == OCMAT_FMT_READABLE); - - // In the compact collinear/SOC layout the header carries the - // per-atom magnetism, which is only known after the Pauli - // block traces; that header is written in the spin branches. - const bool defer_header = readable - && (nspin == 1 || nspin == 2 || nspin == 4); - if (!defer_header) + const int nm = 2 * l + 1; + double electron_count[2] = {0.0, 0.0}; + for (int spin = 0; spin < 2; ++spin) { - if (readable) - { - // Readable header, e.g. " Fe Atom 1 L 2". - // The position is not repeated: it is recorded in - // STRU and constant within one ionic (g) step. - ofs << "\n " << elem_label; - ofs << " Atom " << iat + 1; - ofs << " L " << l << std::endl; - } - else + for (int m = 0; m < nm; ++m) { - // Legacy diagnostic layout. The compact-only restart - // reader intentionally rejects this representation. - ofs << "\n" << elem_label << " " << index_in_type; - ofs << " Atom=" << iat + 1; - ofs << " L=" << l; - ofs << std::setprecision(8) << std::fixed - << std::setw(14) << tau.x - << std::setw(14) << tau.y - << std::setw(14) << tau.z << std::endl; + electron_count[spin] += dftu.occmat().get(iat, l, spin, m, m); } } - if (nspin == 1 || nspin == 2) - { - const int nm = 2 * l + 1; - double sum0[2] = {0.0, 0.0}; + ofs << "\n " << elem_label; + ofs << " Atom " << iat + 1; + ofs << " L " << l; + ofs << " mag " << std::fixed << std::setprecision(8) + << electron_count[0] - electron_count[1] << std::endl; - if (readable) + ofs << std::fixed << std::setprecision(8); + for (int spin = 0; spin < 2; ++spin) + { + ofs << " spin " << spin + 1; + ofs << " nelec " << electron_count[spin] << std::endl; + for (int m0 = 0; m0 < nm; ++m0) { - // Compact snapshot layout, e.g. - // " Fe Atom 1 L 2 mag 0.84215429" - // " spin 1 nelec 1.76153075" - // followed by the matrix rows. - if (diag) - { - for (int is = 0; is < 2; is++) - { - for (int m0 = 0; m0 < nm; m0++) - { - sum0[is] += dftu.occmat().get(iat, l, is, m0, m0); - } - } - } - ofs << "\n " << elem_label; - ofs << " Atom " << iat + 1; - ofs << " L " << l; - if (diag) + for (int m1 = 0; m1 < nm; ++m1) { - ofs << " mag " << std::fixed << std::setprecision(8) - << sum0[0] - sum0[1]; + ofs << " " << std::setw(10) + << dftu.occmat().get(iat, l, spin, m0, m1); } ofs << std::endl; - - ofs << std::fixed << std::setprecision(8); - for (int is = 0; is < 2; is++) - { - ofs << " spin " << is + 1; - if (diag) - { - ofs << " nelec " << sum0[is]; - } - ofs << std::endl; - for (int m0 = 0; m0 < nm; m0++) - { - for (int m1 = 0; m1 < nm; m1++) - { - // One separating space before every - // value; positive values are 10 chars, - // negative ones 11. - ofs << " " << std::setw(10) - << dftu.occmat().get(iat, l, is, m0, m1); - } - ofs << std::endl; - } - } - } - else - { - const std::string trace_label = " Trace (electrons) "; - const std::string matrix_label = " spin= "; - for (int is = 0; is < 2; is++) - { - if (diag) - { - std::vector> A(nm, - std::vector(nm)); - for (int m0 = 0; m0 < nm; m0++) - { - for (int m1 = 0; m1 < nm; m1++) - { - A[m0][m1] = dftu.occmat().get(iat, l, is, m0, m1); - } - } - std::vector eigenvalues = CalculateEigenvalues(A, nm); - sum0[is] = 0.0; - ofs << " Eigenvalues for spin=" << is + 1 << std::endl; - ofs << std::setprecision(8) << std::fixed; - for (int i = 0; i < nm; i++) - { - ofs << std::setw(12) << eigenvalues[i]; - sum0[is] += eigenvalues[i]; - } - ofs << std::endl; - ofs << std::fixed << std::setprecision(8); - ofs << trace_label << sum0[is] << std::endl; - } - ofs << matrix_label << is + 1 << std::endl; - // Matrix elements use fixed-point notation with 8 - // fractional digits (~1e-8 absolute precision). - ofs << std::fixed << std::setprecision(8); - for (int m0 = 0; m0 < nm; m0++) - { - for (int m1 = 0; m1 < nm; m1++) - { - // Legacy keeps its historical - // 12-character column alignment. - ofs << " " << std::setw(11) - << dftu.occmat().get(iat, l, is, m0, m1); - } - ofs << std::endl; - } - } - if (diag) - { - ofs << " mag " << std::fixed << std::setprecision(8) - << sum0[0] - sum0[1] << std::endl; - } } } - else if (nspin == 4) // SOC - { + } + else if (nspin == 4) // SOC + { const int m = 2 * l + 1; std::vector> blocks; extract_soc_pauli_blocks(dftu, iat, l, soc_layout, blocks); - if (readable) - { // Compact snapshot layout aligned with the // nspin == 2 style, e.g. // " Fe Atom 1 L 2 mag mx my mz" @@ -793,22 +662,14 @@ void write_occup_m(const Plus_U_Base& dftu, ofs << "\n " << elem_label; ofs << " Atom " << iat + 1; ofs << " L " << l; - if (diag) - { - ofs << " mag " << std::fixed << std::setprecision(8) - << mag_x << " " << mag_y << " " << mag_z; - } + ofs << " mag " << std::fixed << std::setprecision(8) + << mag_x << " " << mag_y << " " << mag_z; ofs << std::endl; ofs << std::fixed << std::setprecision(8); // spin 1: spin up/up block (Re). - ofs << " spin 1"; - if (diag) - { - ofs << " nelec " << n_up; - } - ofs << std::endl; + ofs << " spin 1 nelec " << n_up << std::endl; for (int m0 = 0; m0 < m; ++m0) { for (int m1 = 0; m1 < m; ++m1) @@ -847,12 +708,7 @@ void write_occup_m(const Plus_U_Base& dftu, } // spin 2: spin down/down block (Re). - ofs << " spin 2"; - if (diag) - { - ofs << " nelec " << n_down; - } - ofs << std::endl; + ofs << " spin 2 nelec " << n_down << std::endl; for (int m0 = 0; m0 < m; ++m0) { for (int m1 = 0; m1 < m; ++m1) @@ -863,130 +719,9 @@ void write_occup_m(const Plus_U_Base& dftu, } ofs << std::endl; } - } - else if (diag) - { - // Occupation numbers of the correlated spinor shell. - std::vector eigenvalues = calculate_spinor_eigenvalues(blocks, m); - ofs << " Eigenvalues in spinor space:" << std::endl; - ofs << std::setprecision(8) << std::fixed; - for (int i = 0; i < 2 * m; ++i) - { - ofs << std::setw(12) << eigenvalues[i]; - if ((i + 1) % m == 0) - { - ofs << std::endl; - } - } - - // Per-atom charge and spin-moment traces of the shell. - double n_tot = 0.0; - double mag_x = 0.0; - double mag_y = 0.0; - double mag_z = 0.0; - for (int mm = 0; mm < m; ++mm) - { - const int k = mm * m + mm; - n_tot += blocks[0][k]; - mag_x += blocks[1][k]; - mag_y += blocks[2][k]; - mag_z += blocks[3][k]; - } - ofs << std::fixed << std::setprecision(8); - ofs << " Trace (electrons) " << n_tot << std::endl; - ofs << " Magnetism for atom " << iat + 1 << " (mx, my, mz):" - << std::setw(12) << mag_x << std::setw(12) << mag_y - << std::setw(12) << mag_z << std::endl; - - // Matrix elements use fixed-point notation with 8 - // fractional digits (~1e-8 absolute precision). - ofs << std::fixed << std::setprecision(8); - - // Hermitian occupation matrix shown spin block by spin block. - ofs << " Occupation matrix, spin up/up (Re):" << std::endl; - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double val = 0.5 * (blocks[0][k] + blocks[3][k]); - ofs << " " << std::setw(10) << val; - } - ofs << std::endl; - } - ofs << " Occupation matrix, spin up/down (Re):" << std::endl; - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double val = 0.5 * blocks[1][k]; - ofs << " " << std::setw(10) << val; - } - ofs << std::endl; - } - ofs << " Occupation matrix, spin up/down (Im):" << std::endl; - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double val = 0.5 * blocks[2][k]; - ofs << " " << std::setw(10) << val; - } - ofs << std::endl; - } - ofs << " Occupation matrix, spin down/down (Re):" << std::endl; - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double val = 0.5 * (blocks[0][k] - blocks[3][k]); - ofs << " " << std::setw(10) << val; - } - ofs << std::endl; - } - } - else - { - // Legacy real 2m x 2m matrix in spin basis; retained - // for diagnostic output only. - for (int m0 = 0; m0 < m; ++m0) - { - for (int ipol0 = 0; ipol0 < npol; ++ipol0) - { - for (int m1 = 0; m1 < m; ++m1) - { - for (int ipol1 = 0; ipol1 < npol; ++ipol1) - { - const int k = m0 * m + m1; - double val = 0.0; - if (ipol0 == 0 && ipol1 == 0) - { - val = 0.5 * (blocks[0][k] + blocks[3][k]); - } - else if (ipol0 == 1 && ipol1 == 1) - { - val = 0.5 * (blocks[0][k] - blocks[3][k]); - } - else - { - val = 0.5 * blocks[1][k]; - } - ofs << " " << std::setw(10) - << std::setprecision(8) - << std::fixed << val; - } - } - ofs << std::endl; - } - } - } - } - } // l - } // I - } // T + } + } + } return; } diff --git a/source/source_pw/module_pwdft/dftu_base_io.h b/source/source_pw/module_pwdft/dftu_base_io.h index babcc76908..5a43b1debf 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.h +++ b/source/source_pw/module_pwdft/dftu_base_io.h @@ -189,19 +189,15 @@ void write_latest_occmat(const Plus_U_Base& dftu, /// Write local occupation matrices to the given stream. /// -/// When diag is true, eigenvalues and magnetism are also printed; otherwise -/// only raw matrix elements. fmt selects between the legacy token layout -/// (used for running log) and the readable snapshot layout of -/// occ_mat.txt and occ_matg{#}.txt. soc_layout tells how the nspin == 4 -/// storage is arranged (PW Pauli blocks or LCAO real spin-basis matrix). +/// Writes the compact snapshot layout accepted by read_occup_m(). The output +/// includes per-atom magnetism, electron counts and the full matrix blocks. +/// soc_layout tells how the nspin == 4 storage is arranged (PW Pauli blocks +/// or LCAO real spin-basis matrix). /// Caller is responsible for opening/closing the stream. void write_occup_m(const Plus_U_Base& dftu, const UnitCell& ucell, std::ostream& ofs, - bool diag, int nspin, - int npol, - OccmatTextFormat fmt, OccmatSocLayout soc_layout); } // namespace DFTU_BASE diff --git a/source/source_pw/module_pwdft/dftu_types.h b/source/source_pw/module_pwdft/dftu_types.h index e1f1bee13d..36ffac42a4 100644 --- a/source/source_pw/module_pwdft/dftu_types.h +++ b/source/source_pw/module_pwdft/dftu_types.h @@ -23,13 +23,6 @@ struct OccmatOutputCfg int dft_plus_u; }; -/// Text format used by write_occup_m(). -enum OccmatTextFormat -{ - OCMAT_FMT_LEGACY, ///< legacy diagnostic layout; not accepted on input - OCMAT_FMT_READABLE ///< compact human-readable layout for snapshot files -}; - /// Storage layout of the real (2l+1)*npol matrix when nspin == 4. enum OccmatSocLayout { diff --git a/source/source_pw/module_pwdft/test/dftu_base_test.cpp b/source/source_pw/module_pwdft/test/dftu_base_test.cpp index e5197b98eb..01a24fbb60 100644 --- a/source/source_pw/module_pwdft/test/dftu_base_test.cpp +++ b/source/source_pw/module_pwdft/test/dftu_base_test.cpp @@ -865,10 +865,7 @@ TEST_F(DFTUBaseTest, CompactWriterReaderCollinearRoundtrip) DFTU_BASE::write_occup_m(source, ucell, output, - true, 2, - 1, - DFTU_BASE::OCMAT_FMT_READABLE, DFTU_BASE::SOC_LAYOUT_PAULI); } @@ -918,10 +915,7 @@ TEST_F(DFTUBaseTest, CompactWriterReaderSocPauliRoundtrip) DFTU_BASE::write_occup_m(source, ucell, output, - true, 4, - 2, - DFTU_BASE::OCMAT_FMT_READABLE, DFTU_BASE::SOC_LAYOUT_PAULI); } From 37f7cabb9784030a44a689b87a82606776c95edc Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Wed, 21 Oct 2026 22:47:27 +0800 Subject: [PATCH 3/6] Split DFT+U compact occupation writer --- source/Makefile.Objects | 1 + .../module_dftu/unittests/CMakeLists.txt | 2 + source/source_pw/module_pwdft/CMakeLists.txt | 1 + .../source_pw/module_pwdft/dftu_base_io.cpp | 343 ------------------ .../module_pwdft/dftu_base_io_write.cpp | 329 +++++++++++++++++ .../module_pwdft/test/CMakeLists.txt | 1 + .../module_pwdft/test/dftu_base_test.cpp | 52 +++ 7 files changed, 386 insertions(+), 343 deletions(-) create mode 100644 source/source_pw/module_pwdft/dftu_base_io_write.cpp diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 191f5e4b78..6d2615eddf 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -958,6 +958,7 @@ OBJS_SRCPW=h_ewald_pw.o\ dftu_base.o\ dftu_base_io.o\ dftu_base_io_read.o\ + dftu_base_io_write.o\ dftu_pw.o\ dftu_pw_tools.o\ yukawa_screening.o\ diff --git a/source/source_lcao/module_dftu/unittests/CMakeLists.txt b/source/source_lcao/module_dftu/unittests/CMakeLists.txt index 6a50a09ffe..fb94f80174 100644 --- a/source/source_lcao/module_dftu/unittests/CMakeLists.txt +++ b/source/source_lcao/module_dftu/unittests/CMakeLists.txt @@ -92,6 +92,7 @@ AddTest( ../../../source_pw/module_pwdft/dftu_base.cpp ../../../source_pw/module_pwdft/dftu_base_io.cpp ../../../source_pw/module_pwdft/dftu_base_io_read.cpp + ../../../source_pw/module_pwdft/dftu_base_io_write.cpp ../../../source_pw/module_pwdft/yukawa_screening.cpp ../../../source_estate/occ_matrix.cpp ../../../source_estate/occ_mixer.cpp @@ -138,6 +139,7 @@ AddTest( ../../../source_pw/module_pwdft/dftu_base.cpp ../../../source_pw/module_pwdft/dftu_base_io.cpp ../../../source_pw/module_pwdft/dftu_base_io_read.cpp + ../../../source_pw/module_pwdft/dftu_base_io_write.cpp ../../../source_pw/module_pwdft/yukawa_screening.cpp ../../../source_estate/occ_matrix.cpp ../../../source_estate/occ_mixer.cpp diff --git a/source/source_pw/module_pwdft/CMakeLists.txt b/source/source_pw/module_pwdft/CMakeLists.txt index df29445f84..b78a019700 100644 --- a/source/source_pw/module_pwdft/CMakeLists.txt +++ b/source/source_pw/module_pwdft/CMakeLists.txt @@ -34,6 +34,7 @@ list(APPEND objects dftu_base.cpp dftu_base_io.cpp dftu_base_io_read.cpp + dftu_base_io_write.cpp dftu_pw.cpp dftu_pw_tools.cpp yukawa_screening.cpp diff --git a/source/source_pw/module_pwdft/dftu_base_io.cpp b/source/source_pw/module_pwdft/dftu_base_io.cpp index db5ffca1fc..38e3bddcb8 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.cpp +++ b/source/source_pw/module_pwdft/dftu_base_io.cpp @@ -4,200 +4,17 @@ #include "source_pw/module_pwdft/dftu_base.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" -#include "source_base/module_external/lapack_connector.h" #include "source_base/parallel_common.h" #include "source_base/parallel_global.h" #include "source_base/timer.h" #include "source_main/version.h" -#include -#include -#include #include #include #include #include #include -namespace -{ - -// 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 file orchestration. -inline void JacobiRotate(std::vector>& A, int p, int q, int n) -{ - if (std::abs(A[p][q]) > 1e-10) - { - double r = (A[q][q] - A[p][p]) / (2.0 * A[p][q]); - double t = 0.0; - if (r >= 0) - { - t = 1.0 / (r + sqrt(1.0 + r * r)); - } - else - { - t = -1.0 / (-r + sqrt(1.0 + r * r)); - } - double c = 1.0 / sqrt(1.0 + t * t); - double s = t * c; - - A[p][p] -= t * A[p][q]; - A[q][q] += t * A[p][q]; - A[p][q] = A[q][p] = 0.0; - - for (int k = 0; k < n; k++) - { - if (k != p && k != q) - { - double Akp = c * A[k][p] - s * A[k][q]; - double Akq = s * A[k][p] + c * A[k][q]; - A[k][p] = A[p][k] = Akp; - A[k][q] = A[q][k] = Akq; - } - } - } -} - -inline std::vector CalculateEigenvalues(std::vector>& A, int n) -{ - std::vector eigenvalues(n); - while (true) - { - int p = 0, q = 1; - for (int i = 0; i < n; i++) - { - for (int j = i + 1; j < n; j++) - { - if (std::abs(A[i][j]) > std::abs(A[p][q])) - { - p = i; - q = j; - } - } - } - - if (std::abs(A[p][q]) < 1e-10) - { - for (int i = 0; i < n; i++) - { - eigenvalues[i] = A[i][i]; - } - break; - } - - JacobiRotate(A, p, q, n); - } - return eigenvalues; -} - -/// @brief Extract the four Pauli-component blocks of one SOC correlated shell. -/// -/// The blocks, each (2l+1)x(2l+1) stored row-major with index m0*m+m1, are -/// the charge channel and the sigma_x/y/z spin channels. The PW storage keeps -/// the four blocks contiguous; the LCAO storage is a real symmetric matrix in -/// spin basis, whose spin-off-diagonal imaginary channel is not stored. -/// -/// @param dftu DFT+U object holding the occupation matrix -/// @param iat atom index -/// @param l angular-momentum channel -/// @param layout storage layout of the nspin == 4 occupation matrix -/// @param blocks [out] four Pauli-component blocks -inline void extract_soc_pauli_blocks(const Plus_U_Base& dftu, - const int iat, - const int l, - const DFTU_BASE::OccmatSocLayout layout, - std::vector>& blocks) -{ - const int m = 2 * l + 1; - const int m2 = m * m; - blocks.assign(4, std::vector(m2, 0.0)); - - const OccupationMatrix& occmat = dftu.occmat(); - if (layout == DFTU_BASE::SOC_LAYOUT_PAULI) - { - const ModuleBase::matrix& occ = occmat.mat(iat, l, 0); - for (int is = 0; is < 4; ++is) - { - for (int k = 0; k < m2; ++k) - { - blocks[is][k] = occ.c[is * m2 + k]; - } - } - } - else // SOC_LAYOUT_SPIN_BASIS_REAL - { - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double suu = occmat.get(iat, l, 0, m0, m1); - const double sud = occmat.get(iat, l, 0, m0, m1 + m); - const double sdu = occmat.get(iat, l, 0, m0 + m, m1); - const double sdd = occmat.get(iat, l, 0, m0 + m, m1 + m); - blocks[0][k] = suu + sdd; - blocks[1][k] = sud + sdu; - // The spin-basis storage drops imaginary parts, so the - // sigma_y channel has no contribution. - blocks[2][k] = 0.0; - blocks[3][k] = suu - sdd; - } - } - } -} - -inline std::vector calculate_spinor_eigenvalues( - const std::vector>& blocks, - const int m) -{ - const int n = 2 * m; - ModuleBase::ComplexMatrix spinor_n(n, n); - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double b0 = blocks[0][k]; - const double b1 = blocks[1][k]; - const double b2 = blocks[2][k]; - const double b3 = blocks[3][k]; - - spinor_n(m0, m1) = std::complex(0.5 * (b0 + b3), 0.0); - spinor_n(m0, m1 + m) = std::complex(0.5 * b1, 0.5 * b2); - spinor_n(m0 + m, m1) = std::complex(0.5 * b1, -0.5 * b2); - spinor_n(m0 + m, m1 + m) = std::complex(0.5 * (b0 - b3), 0.0); - } - } - - std::vector eigenvalues(n, 0.0); - std::vector> work(1, std::complex(0.0, 0.0)); - std::vector rwork(std::max(1, 3 * n - 2), 0.0); - int info = 0; - int lwork = -1; - - LapackConnector::zheev('N', 'U', n, spinor_n, n, eigenvalues.data(), - work.data(), lwork, rwork.data(), &info); - if (info != 0) - { - ModuleBase::WARNING_QUIT("calculate_spinor_eigenvalues", - "LAPACK zheev workspace query failed"); - } - lwork = std::max(1, static_cast(work[0].real())); - work.resize(lwork); - LapackConnector::zheev('N', 'U', n, spinor_n, n, eigenvalues.data(), - work.data(), lwork, rwork.data(), &info); - if (info != 0) - { - ModuleBase::WARNING_QUIT("calculate_spinor_eigenvalues", - "LAPACK zheev diagonalization failed"); - } - return eigenvalues; -} - -} // namespace - - namespace { @@ -567,164 +384,4 @@ void write_latest_occmat(const Plus_U_Base& dftu, } -void write_occup_m(const Plus_U_Base& dftu, - const UnitCell& ucell, - std::ostream& ofs, - int nspin, - OccmatSocLayout soc_layout) -{ - ModuleBase::TITLE("DFTU_BASE", "write_occup_m"); - - if (GlobalV::MY_RANK != 0) - { - return; - } - - for (int T = 0; T < ucell.ntype; T++) - { - if (!dftu.has_l_channel(T)) - { - continue; - } - const int l = dftu.get_l_channel(T); - - for (int I = 0; I < ucell.atoms[T].na; I++) - { - const int iat = ucell.itia2iat(T, I); - const std::string& elem_label = ucell.atoms[T].label; - - if (nspin == 1 || nspin == 2) - { - const int nm = 2 * l + 1; - double electron_count[2] = {0.0, 0.0}; - for (int spin = 0; spin < 2; ++spin) - { - for (int m = 0; m < nm; ++m) - { - electron_count[spin] += dftu.occmat().get(iat, l, spin, m, m); - } - } - - ofs << "\n " << elem_label; - ofs << " Atom " << iat + 1; - ofs << " L " << l; - ofs << " mag " << std::fixed << std::setprecision(8) - << electron_count[0] - electron_count[1] << std::endl; - - ofs << std::fixed << std::setprecision(8); - for (int spin = 0; spin < 2; ++spin) - { - ofs << " spin " << spin + 1; - ofs << " nelec " << electron_count[spin] << std::endl; - for (int m0 = 0; m0 < nm; ++m0) - { - for (int m1 = 0; m1 < nm; ++m1) - { - ofs << " " << std::setw(10) - << dftu.occmat().get(iat, l, spin, m0, m1); - } - ofs << std::endl; - } - } - } - else if (nspin == 4) // SOC - { - const int m = 2 * l + 1; - std::vector> blocks; - extract_soc_pauli_blocks(dftu, iat, l, soc_layout, blocks); - - // Compact snapshot layout aligned with the - // nspin == 2 style, e.g. - // " Fe Atom 1 L 2 mag mx my mz" - // " spin 1 nelec ..." - // followed by the four spin blocks. - double mag_x = 0.0; - double mag_y = 0.0; - double mag_z = 0.0; - for (int mm = 0; mm < m; ++mm) - { - const int k = mm * m + mm; - mag_x += blocks[1][k]; - mag_y += blocks[2][k]; - mag_z += blocks[3][k]; - } - // Spin-resolved electron counts: traces of the - // up/up and down/down blocks. - double n_up = 0.0; - double n_down = 0.0; - for (int mm = 0; mm < m; ++mm) - { - const int k = mm * m + mm; - n_up += 0.5 * (blocks[0][k] + blocks[3][k]); - n_down += 0.5 * (blocks[0][k] - blocks[3][k]); - } - - ofs << "\n " << elem_label; - ofs << " Atom " << iat + 1; - ofs << " L " << l; - ofs << " mag " << std::fixed << std::setprecision(8) - << mag_x << " " << mag_y << " " << mag_z; - ofs << std::endl; - - ofs << std::fixed << std::setprecision(8); - - // spin 1: spin up/up block (Re). - ofs << " spin 1 nelec " << n_up << std::endl; - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double val = 0.5 * (blocks[0][k] + blocks[3][k]); - ofs << " " << std::setw(10) << val; - } - ofs << std::endl; - } - - // Off-diagonal Re block. - ofs << " spin 12 re" << std::endl; - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double val = 0.5 * blocks[1][k]; - ofs << " " << std::setw(10) << val; - } - ofs << std::endl; - } - - // Off-diagonal Im block. - ofs << " spin 12 im" << std::endl; - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double val = 0.5 * blocks[2][k]; - ofs << " " << std::setw(10) << val; - } - ofs << std::endl; - } - - // spin 2: spin down/down block (Re). - ofs << " spin 2 nelec " << n_down << std::endl; - for (int m0 = 0; m0 < m; ++m0) - { - for (int m1 = 0; m1 < m; ++m1) - { - const int k = m0 * m + m1; - const double val = 0.5 * (blocks[0][k] - blocks[3][k]); - ofs << " " << std::setw(10) << val; - } - ofs << std::endl; - } - } - } - } - - return; -} - - } // namespace DFTU_BASE diff --git a/source/source_pw/module_pwdft/dftu_base_io_write.cpp b/source/source_pw/module_pwdft/dftu_base_io_write.cpp new file mode 100644 index 0000000000..8b409f7946 --- /dev/null +++ b/source/source_pw/module_pwdft/dftu_base_io_write.cpp @@ -0,0 +1,329 @@ +#include "source_pw/module_pwdft/dftu_base_io.h" + +#include "source_base/global_function.h" +#include "source_base/module_external/lapack_connector.h" +#include "source_cell/unitcell.h" +#include "source_pw/module_pwdft/dftu_base.h" + +#include +#include +#include +#include +#include +#include + +namespace DFTU_BASE +{ + +struct SocCompactBlocks +{ + std::vector up_up; + std::vector up_down_real; + std::vector up_down_imag; + 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 +{ + +inline void jacobi_rotate(std::vector>& matrix, + const int p, + const int q, + const int size) +{ + if (std::abs(matrix[p][q]) <= 1.0e-10) + { + return; + } + + const double ratio = (matrix[q][q] - matrix[p][p]) / (2.0 * matrix[p][q]); + const double tangent = ratio >= 0.0 + ? 1.0 / (ratio + std::sqrt(1.0 + ratio * ratio)) + : -1.0 / (-ratio + std::sqrt(1.0 + ratio * ratio)); + const double cosine = 1.0 / std::sqrt(1.0 + tangent * tangent); + const double sine = tangent * cosine; + + matrix[p][p] -= tangent * matrix[p][q]; + matrix[q][q] += tangent * matrix[p][q]; + matrix[p][q] = matrix[q][p] = 0.0; + + for (int k = 0; k < size; ++k) + { + if (k == p || k == q) + { + continue; + } + const double matrix_kp = cosine * matrix[k][p] - sine * matrix[k][q]; + const double matrix_kq = sine * matrix[k][p] + cosine * matrix[k][q]; + matrix[k][p] = matrix[p][k] = matrix_kp; + matrix[k][q] = matrix[q][k] = matrix_kq; + } +} + +inline std::vector calculate_eigenvalues(std::vector>& matrix, + const int size) +{ + std::vector eigenvalues(size); + while (true) + { + int p = 0; + int q = 1; + for (int row = 0; row < size; ++row) + { + for (int column = row + 1; column < size; ++column) + { + if (std::abs(matrix[row][column]) > std::abs(matrix[p][q])) + { + p = row; + q = column; + } + } + } + + if (std::abs(matrix[p][q]) < 1.0e-10) + { + for (int i = 0; i < size; ++i) + { + eigenvalues[i] = matrix[i][i]; + } + break; + } + jacobi_rotate(matrix, p, q, size); + } + return eigenvalues; +} + +inline std::vector calculate_spinor_eigenvalues(const SocCompactBlocks& blocks, + const int orbital_size) +{ + const int spinor_size = 2 * orbital_size; + ModuleBase::ComplexMatrix spinor_matrix(spinor_size, spinor_size); + for (int row = 0; row < orbital_size; ++row) + { + for (int column = 0; column < orbital_size; ++column) + { + const int index = row * orbital_size + column; + spinor_matrix(row, column) = blocks.up_up[index]; + spinor_matrix(row, column + orbital_size) + = std::complex(blocks.up_down_real[index], + blocks.up_down_imag[index]); + spinor_matrix(row + orbital_size, column) + = std::complex(blocks.up_down_real[index], + -blocks.up_down_imag[index]); + spinor_matrix(row + orbital_size, column + orbital_size) + = blocks.down_down[index]; + } + } + + std::vector eigenvalues(spinor_size, 0.0); + std::vector> work(1, std::complex(0.0, 0.0)); + std::vector rwork(std::max(1, 3 * spinor_size - 2), 0.0); + int info = 0; + int lwork = -1; + LapackConnector::zheev('N', + 'U', + spinor_size, + spinor_matrix, + spinor_size, + eigenvalues.data(), + work.data(), + lwork, + rwork.data(), + &info); + if (info != 0) + { + ModuleBase::WARNING_QUIT("calculate_spinor_eigenvalues", + "LAPACK zheev workspace query failed"); + } + + lwork = std::max(1, static_cast(work[0].real())); + work.resize(lwork); + LapackConnector::zheev('N', + 'U', + spinor_size, + spinor_matrix, + spinor_size, + eigenvalues.data(), + work.data(), + lwork, + rwork.data(), + &info); + if (info != 0) + { + ModuleBase::WARNING_QUIT("calculate_spinor_eigenvalues", + "LAPACK zheev diagonalization failed"); + } + return eigenvalues; +} + +} // namespace retained_eigenvalue_helpers + +static SocCompactBlocks extract_soc_compact_blocks(const Plus_U_Base& dftu, + const int iat, + const int l, + const OccmatSocLayout layout) +{ + const int size = 2 * l + 1; + const int block_size = size * size; + SocCompactBlocks blocks; + blocks.up_up.resize(block_size, 0.0); + blocks.up_down_real.resize(block_size, 0.0); + blocks.up_down_imag.resize(block_size, 0.0); + blocks.down_down.resize(block_size, 0.0); + + const OccupationMatrix& occmat = dftu.occmat(); + if (layout == SOC_LAYOUT_PAULI) + { + const ModuleBase::matrix& occ = occmat.mat(iat, l, 0); + for (int index = 0; index < block_size; ++index) + { + const double charge = occ.c[0 * block_size + index]; + const double spin_x = occ.c[1 * block_size + index]; + const double spin_y = occ.c[2 * block_size + index]; + const double spin_z = occ.c[3 * block_size + index]; + blocks.up_up[index] = 0.5 * (charge + spin_z); + blocks.up_down_real[index] = 0.5 * spin_x; + blocks.up_down_imag[index] = 0.5 * spin_y; + blocks.down_down[index] = 0.5 * (charge - spin_z); + } + } + else + { + for (int row = 0; row < size; ++row) + { + for (int column = 0; column < size; ++column) + { + const int index = row * size + column; + const double up_up = occmat.get(iat, l, 0, row, column); + const double up_down = occmat.get(iat, l, 0, row, column + size); + const double down_up = occmat.get(iat, l, 0, row + size, column); + const double down_down = occmat.get(iat, l, 0, row + size, column + size); + blocks.up_up[index] = up_up; + blocks.up_down_real[index] = 0.5 * (up_down + down_up); + blocks.down_down[index] = down_down; + } + } + } + return blocks; +} + +static void write_matrix(std::ostream& os, const double* matrix, const int size) +{ + for (int row = 0; row < size; ++row) + { + for (int column = 0; column < size; ++column) + { + os << " " << std::setw(10) << matrix[row * size + column]; + } + os << std::endl; + } +} + +static void write_collinear_atom(const Plus_U_Base& dftu, + const UnitCell& ucell, + std::ostream& os, + const int type, + const int iat, + const int l) +{ + const int size = 2 * l + 1; + double electron_count[2] = {0.0, 0.0}; + for (int spin = 0; spin < 2; ++spin) + { + for (int m = 0; m < size; ++m) + { + electron_count[spin] += dftu.occmat().get(iat, l, spin, m, m); + } + } + + os << "\n " << ucell.atoms[type].label; + os << " Atom " << iat + 1; + os << " L " << l; + os << " mag " << std::fixed << std::setprecision(8) + << electron_count[0] - electron_count[1] << std::endl; + + os << std::fixed << std::setprecision(8); + for (int spin = 0; spin < 2; ++spin) + { + os << " spin " << spin + 1; + os << " nelec " << electron_count[spin] << std::endl; + write_matrix(os, dftu.occmat().mat(iat, l, spin).c, size); + } +} + +static void write_soc_atom(const Plus_U_Base& dftu, + const UnitCell& ucell, + std::ostream& os, + const int type, + const int iat, + const int l, + const OccmatSocLayout layout) +{ + const int size = 2 * l + 1; + const SocCompactBlocks blocks = extract_soc_compact_blocks(dftu, iat, l, layout); + double electron_up = 0.0; + double electron_down = 0.0; + double magnetization_x = 0.0; + double magnetization_y = 0.0; + double magnetization_z = 0.0; + for (int m = 0; m < size; ++m) + { + const int index = m * size + m; + electron_up += blocks.up_up[index]; + electron_down += blocks.down_down[index]; + magnetization_x += 2.0 * blocks.up_down_real[index]; + magnetization_y += 2.0 * blocks.up_down_imag[index]; + magnetization_z += blocks.up_up[index] - blocks.down_down[index]; + } + + os << "\n " << ucell.atoms[type].label; + os << " Atom " << iat + 1; + os << " L " << l; + os << " mag " << std::fixed << std::setprecision(8) + << magnetization_x << " " << magnetization_y << " " << magnetization_z + << std::endl; + os << std::fixed << std::setprecision(8); + + os << " spin 1 nelec " << electron_up << std::endl; + write_matrix(os, blocks.up_up.data(), size); + os << " spin 12 re" << std::endl; + write_matrix(os, blocks.up_down_real.data(), size); + os << " spin 12 im" << std::endl; + write_matrix(os, blocks.up_down_imag.data(), size); + os << " spin 2 nelec " << electron_down << std::endl; + write_matrix(os, blocks.down_down.data(), size); +} + +void write_occup_m(const Plus_U_Base& dftu, + const UnitCell& ucell, + std::ostream& os, + const int nspin, + const OccmatSocLayout soc_layout) +{ + ModuleBase::TITLE("DFTU_BASE", "write_occup_m"); + for (int type = 0; type < ucell.ntype; ++type) + { + if (!dftu.has_l_channel(type)) + { + continue; + } + const int l = dftu.get_l_channel(type); + for (int atom = 0; atom < ucell.atoms[type].na; ++atom) + { + const int iat = ucell.itia2iat(type, atom); + if (nspin == 4) + { + write_soc_atom(dftu, ucell, os, type, iat, l, soc_layout); + } + else + { + write_collinear_atom(dftu, ucell, os, type, iat, l); + } + } + } +} + +} // namespace DFTU_BASE diff --git a/source/source_pw/module_pwdft/test/CMakeLists.txt b/source/source_pw/module_pwdft/test/CMakeLists.txt index bcbb411578..1b0c556dc7 100644 --- a/source/source_pw/module_pwdft/test/CMakeLists.txt +++ b/source/source_pw/module_pwdft/test/CMakeLists.txt @@ -29,6 +29,7 @@ AddTest( ../dftu_base.cpp ../dftu_base_io.cpp ../dftu_base_io_read.cpp + ../dftu_base_io_write.cpp ../dftu_pw_tools.cpp ../yukawa_screening.cpp ../../../source_estate/occ_matrix.cpp diff --git a/source/source_pw/module_pwdft/test/dftu_base_test.cpp b/source/source_pw/module_pwdft/test/dftu_base_test.cpp index 01a24fbb60..36cd6ea112 100644 --- a/source/source_pw/module_pwdft/test/dftu_base_test.cpp +++ b/source/source_pw/module_pwdft/test/dftu_base_test.cpp @@ -938,6 +938,58 @@ TEST_F(DFTUBaseTest, CompactWriterReaderSocPauliRoundtrip) } } +TEST_F(DFTUBaseTest, CompactWriterReaderSocRealSpinBasisRoundtrip) +{ + Plus_U_Base source; + init_dftu_for_spin(source, false, 4, 2, DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); + const int orbital_size = 5; + for (int row = 0; row < orbital_size; ++row) + { + for (int column = 0; column < orbital_size; ++column) + { + const double up_up = 0.1 + 0.01 * row + 0.001 * column; + const double up_down = 0.2 + 0.01 * row + 0.001 * column; + const double down_down = 0.3 + 0.01 * row + 0.001 * column; + source.occmat().set(0, 2, 0, row, column, up_up); + source.occmat().set(0, 2, 0, row, column + orbital_size, up_down); + source.occmat().set(0, 2, 0, row + orbital_size, column, up_down); + source.occmat().set(0, 2, 0, row + orbital_size, column + orbital_size, down_down); + } + } + + const std::string filename = testing::TempDir() + "compact_writer_soc_real.txt"; + { + std::ofstream output(filename.c_str()); + ASSERT_TRUE(output.is_open()); + DFTU_BASE::write_occup_m(source, + ucell, + output, + 4, + DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); + } + + OccupationMatrix restored; + const std::vector l_channel = {2}; + restored.init(ucell, l_channel, 4, 2); + DFTU_BASE::read_occup_m(ucell, + restored, + l_channel, + filename, + 4, + 2, + DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); + + for (int row = 0; row < 2 * orbital_size; ++row) + { + for (int column = 0; column < 2 * orbital_size; ++column) + { + EXPECT_NEAR(restored.get(0, 2, 0, row, column), + source.occmat().get(0, 2, 0, row, column), + 1e-8); + } + } +} + TEST_F(DFTUBaseTest, ReadOccMatCompactHandlesLongComment) { const std::string fn = testing::TempDir() + "compact_long_comment.txt"; From 2ee0fee029d32889add758f452e1fd76c6353141 Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Wed, 21 Oct 2026 22:47:42 +0800 Subject: [PATCH 4/6] Simplify DFT+U occupation output interfaces --- source/source_esolver/esolver_ks_lcao.cpp | 20 +- source/source_esolver/esolver_ks_pw.cpp | 22 +- source/source_lcao/setup_dftu_lcao.cpp | 6 +- source/source_lcao/setup_dftu_lcao.h | 2 - source/source_pw/module_pwdft/dftu_base.cpp | 4 +- .../source_pw/module_pwdft/dftu_base_io.cpp | 246 ++++++------------ source/source_pw/module_pwdft/dftu_base_io.h | 79 +----- .../module_pwdft/dftu_base_io_write.cpp | 2 +- source/source_pw/module_pwdft/dftu_types.h | 9 + .../source_pw/module_pwdft/setup_dftu_pw.cpp | 18 +- source/source_pw/module_pwdft/setup_dftu_pw.h | 3 +- .../module_pwdft/test/dftu_base_test.cpp | 57 +--- 12 files changed, 135 insertions(+), 333 deletions(-) diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index e6c047a40c..d6137e4d8b 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -677,7 +677,6 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& // Read globalv values once into locals so the repeated DFT+U call sites // below do not each re-enter the global dependency surface. const std::string& global_out_dir = PARAM.globalv.global_out_dir; - const int npol = PARAM.globalv.npol; const bool gamma_only_local = PARAM.globalv.gamma_only_local; finish_dftu_lcao(conv_esolver, this->inp_->dft_plus_u, @@ -689,7 +688,6 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& hamilt_lcao, global_out_dir, this->inp_->nspin, - npol, gamma_only_local, istep, iter, @@ -719,15 +717,14 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& // overwrite occ_mat.txt with the latest occupation matrix and the // actual charge-density residual; ESolver_KS::iter_finish() above has // computed drho for this electronic step. + const DFTU_BASE::OccmatStepInfo occmat_step{istep, + iter, + this->scf_thr, + this->drho}; DFTU_BASE::write_latest_occmat(*this->dftu_, ucell, global_out_dir, - this->inp_->nspin, - npol, - istep, - iter, - this->scf_thr, - this->drho, + occmat_step, occmat_cfg, DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); @@ -737,14 +734,9 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& DFTU_BASE::append_ion_step_snapshot(*this->dftu_, ucell, global_out_dir, - this->inp_->nspin, - npol, - istep, - iter, + occmat_step, conv_esolver, true, - this->scf_thr, - this->drho, occmat_cfg, DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index b77f95c9b0..1a10a37c31 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -249,8 +249,7 @@ void ESolver_KS_PW::iter_init(UnitCell& ucell, const int istep, const init_occmat_cfg, this->kv.isk.data(), this->inp_->nspin, - this->inp_->kpar, - PARAM.globalv.npol); + this->inp_->kpar); // mohan add 2025-11: push DFT+U energy from Plus_U instance to ElecState if (this->inp_->dft_plus_u) @@ -380,22 +379,20 @@ void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int this->inp_->scf_nmax, this->inp_->out_occ_mat, this->inp_->dft_plus_u}; + const DFTU_BASE::OccmatStepInfo occmat_step{istep, + iter, + this->scf_thr, + this->drho}; // Read globalv values once into locals so the repeated DFT+U call sites // below do not each re-enter the global dependency surface. const std::string& global_out_dir = PARAM.globalv.global_out_dir; - const int npol = PARAM.globalv.npol; const bool latest_ready = (iter > 1 || istep > 0); if (latest_ready) { DFTU_BASE::write_latest_occmat(*this->dftu_, ucell, global_out_dir, - this->inp_->nspin, - npol, - istep, - iter, - this->scf_thr, - this->drho, + occmat_step, occmat_cfg, DFTU_BASE::SOC_LAYOUT_PAULI); } @@ -408,14 +405,9 @@ void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int DFTU_BASE::append_ion_step_snapshot(*this->dftu_, ucell, global_out_dir, - this->inp_->nspin, - npol, - istep, - iter, + occmat_step, conv_esolver, occmat_ready, - this->scf_thr, - this->drho, occmat_cfg, DFTU_BASE::SOC_LAYOUT_PAULI); diff --git a/source/source_lcao/setup_dftu_lcao.cpp b/source/source_lcao/setup_dftu_lcao.cpp index 77fe5e0f5d..5e3f7f33c3 100644 --- a/source/source_lcao/setup_dftu_lcao.cpp +++ b/source/source_lcao/setup_dftu_lcao.cpp @@ -49,7 +49,6 @@ void finish_dftu_lcao(const bool conv_esolver, void* hamilt_lcao, const std::string& global_out_dir, int nspin, - int npol, const bool gamma_only_local, int istep, int iter, @@ -82,8 +81,7 @@ void finish_dftu_lcao(const bool conv_esolver, DFTU_LCAO::cal_energy_correction(*dftu_ptr, ucell, PARAM.inp.nspin); } } - DFTU_BASE::output(*dftu_ptr, ucell, global_out_dir, nspin, npol, - istep, iter, occmat_cfg, DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); + DFTU_BASE::prepare_ion_step_output(global_out_dir, istep, iter, occmat_cfg); /// use the converged occupation matrix for next MD/Relax SCF calculation if (conv_esolver) @@ -103,7 +101,6 @@ template void finish_dftu_lcao(const bool conv_esolver, void* hamilt_lcao, const std::string& global_out_dir, int nspin, - int npol, const bool gamma_only_local, int istep, int iter, @@ -119,7 +116,6 @@ template void finish_dftu_lcao>(const bool conv_esolver, void* hamilt_lcao, const std::string& global_out_dir, int nspin, - int npol, const bool gamma_only_local, int istep, int iter, diff --git a/source/source_lcao/setup_dftu_lcao.h b/source/source_lcao/setup_dftu_lcao.h index ddd1a79bcd..cd85efa00b 100644 --- a/source/source_lcao/setup_dftu_lcao.h +++ b/source/source_lcao/setup_dftu_lcao.h @@ -56,7 +56,6 @@ void init_dftu_lcao(int dft_plus_u, * @param hamilt_lcao Hamiltonian LCAO object * @param global_out_dir Output directory for occ_mat.txt and occ_matg{#}.txt * @param nspin Number of spin channels (1, 2, or 4) - * @param npol Number of polarizations * @param gamma_only_local Whether only the Gamma point is used for LCAO * @param istep Ionic-step index, starting from 0 * @param iter Electronic-iteration index, starting from 1 @@ -74,7 +73,6 @@ void finish_dftu_lcao(const bool conv_esolver, void* hamilt_lcao, const std::string& global_out_dir, int nspin, - int npol, const bool gamma_only_local, int istep, int iter, diff --git a/source/source_pw/module_pwdft/dftu_base.cpp b/source/source_pw/module_pwdft/dftu_base.cpp index 0bb70b999a..35e33c4581 100644 --- a/source/source_pw/module_pwdft/dftu_base.cpp +++ b/source/source_pw/module_pwdft/dftu_base.cpp @@ -184,7 +184,7 @@ void Plus_U_Base::init_base(UnitCell& cell, DFTU_BASE::read_occup_m(cell, this->occmat_, this->l_channel, readin_fn, nspin, npol, soc_layout); #ifdef __MPI - DFTU_BASE::local_occup_bcast(cell, this->occmat_, this->l_channel, nspin, npol); + DFTU_BASE::local_occup_bcast(cell, this->occmat_, this->l_channel); #endif this->set_occmat_ready(); @@ -216,7 +216,7 @@ void Plus_U_Base::init_base(UnitCell& cell, DFTU_BASE::read_occup_m(cell, this->occmat_, this->l_channel, readin_fn, nspin, npol, soc_layout); #ifdef __MPI - DFTU_BASE::local_occup_bcast(cell, this->occmat_, this->l_channel, nspin, npol); + DFTU_BASE::local_occup_bcast(cell, this->occmat_, this->l_channel); #endif this->set_occmat_ready(); } diff --git a/source/source_pw/module_pwdft/dftu_base_io.cpp b/source/source_pw/module_pwdft/dftu_base_io.cpp index 38e3bddcb8..fa8023dfaf 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.cpp +++ b/source/source_pw/module_pwdft/dftu_base_io.cpp @@ -15,13 +15,13 @@ #include #include -namespace +namespace DFTU_BASE { // Write the STRU-style provenance header shared by occ_mat.txt and the // per-ionic-step snapshot files: ABACUS version, local timestamp and the // 1-based relaxation step this file belongs to. -void write_provenance_header(std::ostream& os, const int istep) +static void write_provenance_header(std::ostream& os, const int istep) { std::time_t now = std::time(nullptr); std::tm tm_local; @@ -30,14 +30,28 @@ void write_provenance_header(std::ostream& os, const int istep) std::strftime(time_buf, sizeof(time_buf), "%Y-%m-%d %H:%M:%S", &tm_local); os << "# ABACUS version: " << VERSION << "\n"; os << "# Written at " << time_buf << "\n"; - os << "# RELAX STEP " << istep + 1 << "\n"; + os << "# RELAX STEP " << istep + 1 << "\n\n"; } -} // namespace - - -namespace DFTU_BASE +static void validate_output_configuration(const char* context, + const int istep, + const int iter, + const OccmatOutputCfg& cfg) { + if (istep < 0) + { + ModuleBase::WARNING_QUIT(context, "istep must be >= 0"); + } + if (iter < 1) + { + ModuleBase::WARNING_QUIT(context, "iter must be >= 1"); + } + if (cfg.out_freq_ion < 0 || cfg.out_freq_elec < 1 || cfg.scf_nmax < 1) + { + ModuleBase::WARNING_QUIT(context, + "invalid occupation-matrix output frequency configuration"); + } +} bool is_ion_step_output_step(const int istep, const OccmatOutputCfg& cfg) { @@ -87,11 +101,10 @@ std::string find_first_existing_file(const std::string& dir, /// element. void local_occup_bcast(const UnitCell& ucell, OccupationMatrix& occ, - const std::vector& l_channel, - int nspin, - int npol) + const std::vector& l_channel) { ModuleBase::TITLE("DFTU_BASE", "local_occup_bcast"); + const int nspin = occ.nspin(); for (int T = 0; T < ucell.ntype; T++) { @@ -99,43 +112,31 @@ void local_occup_bcast(const UnitCell& ucell, { continue; } - + const int l = l_channel[T]; for (int I = 0; I < ucell.atoms[T].na; I++) { const int iat = ucell.itia2iat(T, I); - const int L = l_channel[T]; - - for (int l = 0; l <= ucell.atoms[T].nwl; l++) + if (nspin == 1 || nspin == 2) { - if (l != l_channel[T]) - { - continue; - } - - if (nspin == 1 || nspin == 2) + for (int spin = 0; spin < 2; spin++) { - for (int spin = 0; spin < 2; spin++) - { - Parallel_Common::bcast_double(occ.mat(iat, l, spin).c, - occ.mat(iat, l, spin).nr * occ.mat(iat, l, spin).nc); - } - } - else if (nspin == 4) // SOC - { - Parallel_Common::bcast_double(occ.mat(iat, l, 0).c, - occ.mat(iat, l, 0).nr * occ.mat(iat, l, 0).nc); + Parallel_Common::bcast_double(occ.mat(iat, l, spin).c, + occ.mat(iat, l, spin).nr * occ.mat(iat, l, spin).nc); } } + else if (nspin == 4) // SOC + { + Parallel_Common::bcast_double(occ.mat(iat, l, 0).c, + occ.mat(iat, l, 0).nr * occ.mat(iat, l, 0).nc); + } } } - return; } #endif -void prepare_ion_step_file(const std::string& global_out_dir, - const int istep, - const OccmatOutputCfg& cfg) +static void prepare_ion_step_file(const std::string& global_out_dir, + const int istep) { // Only the root process creates the per-ionic-step file. if (GlobalV::MY_RANK != 0) @@ -155,59 +156,23 @@ void prepare_ion_step_file(const std::string& global_out_dir, } // Provenance header shared with occ_mat.txt. write_provenance_header(ofs_ion_step, istep); - ofs_ion_step.close(); - return; } - -void output(const Plus_U_Base& dftu, - const UnitCell& ucell, - const std::string& global_out_dir, - int nspin, - int npol, - int istep, - int iter, - const OccmatOutputCfg& cfg, - OccmatSocLayout soc_layout) +void prepare_ion_step_output(const std::string& global_out_dir, + const int istep, + const int iter, + const OccmatOutputCfg& cfg) { - ModuleBase::TITLE("DFTU_BASE", "output"); - - if (istep < 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::output", "istep must be >= 0"); - } - if (iter < 1) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::output", "iter must be >= 1"); - } - if (cfg.out_freq_ion < 0 || cfg.out_freq_elec < 1 || cfg.scf_nmax < 1) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::output", - "invalid occupation-matrix output frequency configuration"); - } - - // The per-type U values and the per-atom occupation matrices are not - // dumped into the running log here: for large cells (thousands of - // atoms) this floods the log at every electronic step. The dedicated - // occupation-matrix files (occ_mat.txt, occ_matg{#}.txt) gated by - // out_occ_mat carry the same information instead. - - // occ_mat.txt is not written here: its section records drho, which - // is only available after the electronic solve. write_latest_occmat() - // overwrites it at the iter_finish stage. - - // At the first electronic step of an output ionic step, (re)create the - // per-ionic-step file and write its header. The LCAO path enters here - // directly; the PW path additionally prepares the file from - // iter_init_dftu_pw(), because it skips output() at istep 0 / iter 1. - // Electronic-step sections are appended later by append_ion_step_snapshot(). + ModuleBase::TITLE("DFTU_BASE", "prepare_ion_step_output"); + validate_output_configuration("DFTU_BASE::prepare_ion_step_output", + istep, + iter, + cfg); const bool ion_step_output = is_ion_step_output_step(istep, cfg); if (cfg.out_occ_mat && ion_step_output && iter == 1) { - prepare_ion_step_file(global_out_dir, istep, cfg); + prepare_ion_step_file(global_out_dir, istep); } - - return; } // Write one electronic-step section shared by occ_mat.txt and the @@ -215,20 +180,16 @@ void output(const Plus_U_Base& dftu, // charge-density convergence threshold, the actual charge-density // residual and the full occupation matrices. When occmat_ready is false, // the matrix body is an "N/A" placeholder. -void write_snapshot_section(std::ostream& os, - const Plus_U_Base& dftu, - const UnitCell& ucell, - const int iter, - const double scf_thr, - const double drho, - const bool occmat_ready, - const int nspin, - const int npol, - const OccmatSocLayout soc_layout) +static void write_snapshot_section(std::ostream& os, + const Plus_U_Base& dftu, + const UnitCell& ucell, + const OccmatStepInfo& step, + const bool occmat_ready, + const OccmatSocLayout soc_layout) { - os << "# Electronic step " << iter << "\n"; - os << "# scf_thr " << std::scientific << std::setprecision(8) << scf_thr << "\n"; - os << "# drho " << std::scientific << std::setprecision(8) << drho << "\n"; + os << "# Electronic step " << step.iter << "\n"; + os << "# scf_thr " << std::scientific << std::setprecision(8) << step.scf_thr << "\n"; + os << "# drho " << std::scientific << std::setprecision(8) << step.drho << "\n"; // The PW path has no occupation matrix at istep 0 / iter 1 unless it was // loaded from file; record an explicit placeholder instead of a silent @@ -236,27 +197,23 @@ void write_snapshot_section(std::ostream& os, // intentionally rejected by the strict restart reader. if (!occmat_ready) { - os << "\n N/A\n"; + os << "\n N/A\n\n"; return; } // Full occupation matrices together with per-atom magnetism. - write_occup_m(dftu, ucell, os, nspin, soc_layout); + write_occup_m(dftu, ucell, os, soc_layout); + os << "\n"; } void append_ion_step_snapshot(const Plus_U_Base& dftu, const UnitCell& ucell, const std::string& global_out_dir, - int nspin, - int npol, - int istep, - int iter, - bool conv_esolver, - bool occmat_ready, - double scf_thr, - double drho, + const OccmatStepInfo& step, + const bool conv_esolver, + const bool occmat_ready, const OccmatOutputCfg& cfg, - OccmatSocLayout soc_layout) + const OccmatSocLayout soc_layout) { ModuleBase::TITLE("DFTU_BASE", "append_ion_step_snapshot"); @@ -264,7 +221,7 @@ void append_ion_step_snapshot(const Plus_U_Base& dftu, // (Plus_U_Base::l_channel stays empty); writing here would dereference // a null l_channel.data() in write_occup_m. Match the documented contract // that out_occ_mat only takes effect for DFT+U calculations. Check this - // before validating iter/istep/nspin: when DFT+U is off this function is + // before validating the step: when DFT+U is off this function is // a no-op, so an invalid iter (e.g. from an upstream bug) must not crash // here. Mirrors the guard-first pattern in write_latest_occmat. if (!cfg.out_occ_mat || cfg.dft_plus_u <= 0) @@ -272,51 +229,19 @@ void append_ion_step_snapshot(const Plus_U_Base& dftu, return; } - if (nspin != 1 && nspin != 2 && nspin != 4) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::append_ion_step_snapshot", "nspin must be 1, 2 or 4"); - } - if (istep < 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::append_ion_step_snapshot", "istep must be >= 0"); - } - if (iter < 1) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::append_ion_step_snapshot", "iter must be >= 1"); - } - // out_freq_ion == 0 is the valid default (no numbered file) and is - // handled by the gates below; only negative values are an error. - if (cfg.out_freq_ion < 0 || cfg.out_freq_elec < 1 || cfg.scf_nmax < 1) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::append_ion_step_snapshot", - "invalid occupation-matrix output frequency configuration"); - } + validate_output_configuration("DFTU_BASE::append_ion_step_snapshot", + step.istep, + step.iter, + cfg); - const bool ion_step_output = is_ion_step_output_step(istep, cfg); - const bool elec_trigger = is_elec_snapshot_trigger(iter, conv_esolver, cfg); + const bool ion_step_output = is_ion_step_output_step(step.istep, cfg); + const bool elec_trigger = is_elec_snapshot_trigger(step.iter, conv_esolver, cfg); if (!ion_step_output || !elec_trigger || GlobalV::MY_RANK != 0) { return; } - const std::string ion_step_fn = gen_ion_step_occ_mat_filename(global_out_dir, istep); - - // The provenance header written by prepare_ion_step_file() must not be - // counted as a recorded section: scan the file for an existing section - // marker instead of testing the stream position. - bool section_exists = false; - { - std::ifstream checker(ion_step_fn.c_str()); - std::string line; - while (std::getline(checker, line)) - { - if (line.rfind("# Electronic step", 0) == 0) - { - section_exists = true; - break; - } - } - } + const std::string ion_step_fn = gen_ion_step_occ_mat_filename(global_out_dir, step.istep); std::ofstream ofs_ion_step; ofs_ion_step.open(ion_step_fn, std::ios::app); @@ -326,34 +251,25 @@ void append_ion_step_snapshot(const Plus_U_Base& dftu, "Can't open per-ionic-step occupation-matrix file"); } - // Separate consecutive sections by a blank line; the first section - // directly follows the provenance header. - if (section_exists) + // A direct caller may not have prepared the file at iter 1. Add the + // provenance header when append mode opens an empty file. + if (ofs_ion_step.tellp() == std::streampos(0)) { - ofs_ion_step << "\n"; + write_provenance_header(ofs_ion_step, step.istep); } - write_snapshot_section(ofs_ion_step, dftu, ucell, iter, scf_thr, drho, - occmat_ready, nspin, npol, soc_layout); - - ofs_ion_step.close(); + write_snapshot_section(ofs_ion_step, dftu, ucell, step, occmat_ready, soc_layout); - GlobalV::ofs_running << " Append local occupation matrices of electronic step " << iter + GlobalV::ofs_running << " Append local occupation matrices of electronic step " << step.iter << " to file: " << ion_step_fn << std::endl; - return; } void write_latest_occmat(const Plus_U_Base& dftu, const UnitCell& ucell, const std::string& global_out_dir, - int nspin, - int npol, - int istep, - int iter, - double scf_thr, - double drho, + const OccmatStepInfo& step, const OccmatOutputCfg& cfg, - OccmatSocLayout soc_layout) + const OccmatSocLayout soc_layout) { ModuleBase::TITLE("DFTU_BASE", "write_latest_occmat"); @@ -372,15 +288,11 @@ void write_latest_occmat(const Plus_U_Base& dftu, // occ_mat.txt is a single-section snapshot file: the provenance // header shared with the g files followed by the same section writer. - write_provenance_header(ofdftu, istep); + write_provenance_header(ofdftu, step.istep); // Callers of write_latest_occmat() guarantee the matrix exists: the PW // esolver gates on latest_ready, the LCAO esolver computes it every // electronic iteration. - write_snapshot_section(ofdftu, dftu, ucell, iter, scf_thr, drho, - true, nspin, npol, soc_layout); - - ofdftu.close(); - return; + write_snapshot_section(ofdftu, dftu, ucell, step, true, 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 5a43b1debf..4750cf0fbe 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.h +++ b/source/source_pw/module_pwdft/dftu_base_io.h @@ -56,36 +56,24 @@ std::string find_first_existing_file(const std::string& dir, /// of the current electronic step, followed by the occupation matrices and /// per-atom magnetism. When occmat_ready is false, an "N/A" placeholder is /// recorded instead of the matrix body (the PW path has no matrix at -/// istep 0 / iter 1 unless it was loaded from file). The caller is -/// responsible for truncating the file at the first electronic step of the -/// ionic step (see prepare_ion_step_file()) and for invoking this function -/// only after drho of the current electronic step is computed. +/// istep 0 / iter 1 unless it was loaded from file). The caller invokes this +/// function only after drho of the current electronic step is computed. /// /// @param dftu DFT+U object holding the occupation matrices /// @param ucell unit cell /// @param global_out_dir output directory (including the trailing separator) -/// @param nspin number of spin components (1, 2 or 4) -/// @param npol number of polarizations -/// @param istep ionic-step index, starting from 0 -/// @param iter electronic-iteration index, starting from 1 +/// @param step electronic-step metadata and convergence values /// @param conv_esolver whether the electronic SCF is converged at this step /// @param occmat_ready whether the occupation matrix of this step exists; /// false records an "N/A" placeholder instead of the matrix body -/// @param scf_thr configured charge-density convergence threshold -/// @param drho actual charge-density residual of the current electronic step /// @param cfg frequency configuration /// @param soc_layout storage layout of the nspin == 4 occupation matrix void append_ion_step_snapshot(const Plus_U_Base& dftu, const UnitCell& ucell, const std::string& global_out_dir, - int nspin, - int npol, - int istep, - int iter, + const OccmatStepInfo& step, bool conv_esolver, bool occmat_ready, - double scf_thr, - double drho, const OccmatOutputCfg& cfg, OccmatSocLayout soc_layout); @@ -113,46 +101,18 @@ void read_occup_m(const UnitCell& ucell, /// Implemented in dftu_base_io.cpp (only available in MPI builds). void local_occup_bcast(const UnitCell& ucell, OccupationMatrix& occ, - const std::vector& l_channel, - int nspin, - int npol); + const std::vector& l_channel); /// Create (or truncate) the per-ionic-step file occ_matg{istep+1}.txt and /// write its header (rank 0 only). /// -/// Must be called once at the first electronic step (iter == 1) of an output -/// ionic step, before any append_ion_step_snapshot() section. Truncating here -/// also guarantees that a rerun in the same output directory can not append -/// snapshots of a previous calculation. -void prepare_ion_step_file(const std::string& global_out_dir, - const int istep, - const OccmatOutputCfg& cfg); - -/// Output DFT+U information (Hubbard U/J, local occupation matrices) to the -/// running log. -/// -/// When cfg.out_occ_mat is true, cfg.out_freq_ion is positive and istep is an -/// output ionic step, the per-ionic-step file occ_matg{istep+1}.txt is -/// created: at the first electronic step (iter == 1) it is truncated and -/// initialized with a provenance header, then append_ion_step_snapshot() -/// appends one section per recorded electronic step. -/// -/// Note: occ_mat.txt is NOT written here. It records the actual -/// charge-density residual drho, which is only known after the electronic -/// solve, so write_latest_occmat() writes it from the iter_finish stage. -/// -/// Extracted from Plus_U_Base::output as a free function so that IO logic is -/// decoupled from the Plus_U_Base class. The function only reads the -/// Plus_U_Base state via public accessors; no friend declaration needed. -void output(const Plus_U_Base& dftu, - const UnitCell& ucell, - const std::string& global_out_dir, - int nspin, - int npol, - int istep, - int iter, - const OccmatOutputCfg& cfg, - OccmatSocLayout soc_layout); +/// At iter == 1 of an enabled output ionic step, truncates the file and writes +/// its provenance header. Other calls are no-ops. Truncating here guarantees +/// that a rerun can not append snapshots from a previous calculation. +void prepare_ion_step_output(const std::string& global_out_dir, + int istep, + int iter, + const OccmatOutputCfg& cfg); /// Overwrite occ_mat.txt with the occupation matrix of the current /// electronic step (rank 0 only). @@ -167,23 +127,13 @@ void output(const Plus_U_Base& dftu, /// @param dftu DFT+U object holding the occupation matrices /// @param ucell unit cell /// @param global_out_dir output directory (including the trailing separator) -/// @param nspin number of spin components (1, 2 or 4) -/// @param npol number of polarizations -/// @param istep ionic-step index, starting from 0 -/// @param iter electronic-iteration index, starting from 1 -/// @param scf_thr configured charge-density convergence threshold -/// @param drho actual charge-density residual of the current electronic step +/// @param step electronic-step metadata and convergence values /// @param cfg output configuration; only cfg.out_occ_mat is consulted here /// @param soc_layout storage layout of the nspin == 4 occupation matrix void write_latest_occmat(const Plus_U_Base& dftu, const UnitCell& ucell, const std::string& global_out_dir, - int nspin, - int npol, - int istep, - int iter, - double scf_thr, - double drho, + const OccmatStepInfo& step, const OccmatOutputCfg& cfg, OccmatSocLayout soc_layout); @@ -197,7 +147,6 @@ void write_latest_occmat(const Plus_U_Base& dftu, void write_occup_m(const Plus_U_Base& dftu, const UnitCell& ucell, std::ostream& ofs, - int nspin, OccmatSocLayout soc_layout); } // namespace DFTU_BASE diff --git a/source/source_pw/module_pwdft/dftu_base_io_write.cpp b/source/source_pw/module_pwdft/dftu_base_io_write.cpp index 8b409f7946..d10fe5038c 100644 --- a/source/source_pw/module_pwdft/dftu_base_io_write.cpp +++ b/source/source_pw/module_pwdft/dftu_base_io_write.cpp @@ -300,10 +300,10 @@ static void write_soc_atom(const Plus_U_Base& dftu, void write_occup_m(const Plus_U_Base& dftu, const UnitCell& ucell, std::ostream& os, - const int nspin, const OccmatSocLayout soc_layout) { ModuleBase::TITLE("DFTU_BASE", "write_occup_m"); + const int nspin = dftu.occmat().nspin(); for (int type = 0; type < ucell.ntype; ++type) { if (!dftu.has_l_channel(type)) diff --git a/source/source_pw/module_pwdft/dftu_types.h b/source/source_pw/module_pwdft/dftu_types.h index 36ffac42a4..a2e35b9bd8 100644 --- a/source/source_pw/module_pwdft/dftu_types.h +++ b/source/source_pw/module_pwdft/dftu_types.h @@ -23,6 +23,15 @@ struct OccmatOutputCfg int dft_plus_u; }; +/// Electronic-step metadata recorded in an occupation-matrix snapshot. +struct OccmatStepInfo +{ + int istep; + int iter; + double scf_thr; + double drho; +}; + /// Storage layout of the real (2l+1)*npol matrix when nspin == 4. enum OccmatSocLayout { diff --git a/source/source_pw/module_pwdft/setup_dftu_pw.cpp b/source/source_pw/module_pwdft/setup_dftu_pw.cpp index 6e0332eac8..d3be1285fa 100644 --- a/source/source_pw/module_pwdft/setup_dftu_pw.cpp +++ b/source/source_pw/module_pwdft/setup_dftu_pw.cpp @@ -18,24 +18,15 @@ void iter_init_dftu_pw(const int iter, const OccmatOutputCfg& occmat_cfg, const int* isk, const int nspin, - const int kpar, - const int npol) + const int kpar) { if (!p_chgmix || occmat_cfg.dft_plus_u <= 0) { return; } - // Prepare the per-ionic-step file at the first electronic step, before - // the first occupation matrix exists. output() is skipped below at - // istep 0 / iter 1; without this call the g1 file would miss its header - // and could inherit stale sections from a previous run. - if (iter == 1 && occmat_cfg.out_occ_mat - && DFTU_BASE::is_ion_step_output_step(istep, occmat_cfg)) - { - DFTU_BASE::prepare_ion_step_file(global_out_dir, - istep, occmat_cfg); - } + // Prepare the per-ionic-step file before the first snapshot is appended. + DFTU_BASE::prepare_ion_step_output(global_out_dir, istep, iter, occmat_cfg); const bool first_scf_step = (iter == 1 && istep == 0); @@ -82,9 +73,6 @@ void iter_init_dftu_pw(const int iter, return; } - DFTU_BASE::output(dftu, ucell, global_out_dir, - nspin, npol, istep, iter, occmat_cfg, - DFTU_BASE::SOC_LAYOUT_PAULI); } } diff --git a/source/source_pw/module_pwdft/setup_dftu_pw.h b/source/source_pw/module_pwdft/setup_dftu_pw.h index 05c448e216..0b242ce852 100644 --- a/source/source_pw/module_pwdft/setup_dftu_pw.h +++ b/source/source_pw/module_pwdft/setup_dftu_pw.h @@ -26,8 +26,7 @@ void iter_init_dftu_pw(const int iter, const OccmatOutputCfg& occmat_cfg, const int* isk, const int nspin, - const int kpar, - const int npol); + const int kpar); } diff --git a/source/source_pw/module_pwdft/test/dftu_base_test.cpp b/source/source_pw/module_pwdft/test/dftu_base_test.cpp index 36cd6ea112..decf278af8 100644 --- a/source/source_pw/module_pwdft/test/dftu_base_test.cpp +++ b/source/source_pw/module_pwdft/test/dftu_base_test.cpp @@ -462,20 +462,17 @@ TEST_F(DFTUBaseTest, AppendSnapshotNAPlaceholderAndReady) const DFTU_BASE::OccmatOutputCfg cfg = {1, 1, 5, true, 1}; const std::string out_dir = testing::TempDir(); + const DFTU_BASE::OccmatStepInfo first_step = {0, 1, 1e-6, 0.5}; + const DFTU_BASE::OccmatStepInfo second_step = {1, 2, 1e-6, 0.5}; // fresh run: no occupation matrix has been computed or loaded ASSERT_FALSE(dftu.is_occmat_ready()); DFTU_BASE::append_ion_step_snapshot(dftu, ucell, out_dir, - 2, // nspin - 1, // npol - 0, // istep -> occ_matg1.txt - 1, // iter + first_step, false, false, // occmat_ready - 1e-6, - 0.5, cfg, DFTU_BASE::SOC_LAYOUT_PAULI); @@ -496,14 +493,9 @@ TEST_F(DFTUBaseTest, AppendSnapshotNAPlaceholderAndReady) DFTU_BASE::append_ion_step_snapshot(dftu, ucell, out_dir, - 2, // nspin - 1, // npol - 1, // istep -> occ_matg2.txt - 2, // iter + second_step, false, true, // occmat_ready - 1e-6, - 0.5, cfg, DFTU_BASE::SOC_LAYOUT_PAULI); @@ -532,6 +524,7 @@ TEST_F(DFTUBaseTest, OccMatSwitchDisabledWritesNothing) const DFTU_BASE::OccmatOutputCfg cfg = {1, 1, 5, false, 1}; const std::string out_dir = testing::TempDir(); + const DFTU_BASE::OccmatStepInfo step = {0, 1, 1e-6, 0.5}; // Remove stale files from previous test binary runs so the // "file not written" assertions are meaningful. @@ -541,14 +534,9 @@ TEST_F(DFTUBaseTest, OccMatSwitchDisabledWritesNothing) DFTU_BASE::append_ion_step_snapshot(dftu, ucell, out_dir, - 2, // nspin - 1, // npol - 0, // istep (an output ionic step) - 1, // iter + step, false, true, // occmat_ready - 1e-6, - 0.5, cfg, DFTU_BASE::SOC_LAYOUT_PAULI); @@ -559,12 +547,7 @@ TEST_F(DFTUBaseTest, OccMatSwitchDisabledWritesNothing) DFTU_BASE::write_latest_occmat(dftu, ucell, out_dir, - 2, // nspin - 1, // npol - 0, // istep - 2, // iter - 1e-6, - 0.5, + step, cfg, DFTU_BASE::SOC_LAYOUT_PAULI); @@ -589,18 +572,14 @@ TEST_F(DFTUBaseTest, DftPlusUDisabledWritesNothing) const DFTU_BASE::OccmatOutputCfg cfg = {1, 1, 5, true, 0}; const std::string out_dir = testing::TempDir() + "guard_test/"; + const DFTU_BASE::OccmatStepInfo step = {0, 1, 1e-6, 0.5}; DFTU_BASE::append_ion_step_snapshot(dftu, ucell, out_dir, - 2, // nspin - 1, // npol - 0, // istep (an output ionic step) - 1, // iter + step, false, true, // occmat_ready - 1e-6, - 0.5, cfg, DFTU_BASE::SOC_LAYOUT_PAULI); @@ -611,12 +590,7 @@ TEST_F(DFTUBaseTest, DftPlusUDisabledWritesNothing) DFTU_BASE::write_latest_occmat(dftu, ucell, out_dir, - 2, // nspin - 1, // npol - 0, // istep - 2, // iter - 1e-6, - 0.5, + step, cfg, DFTU_BASE::SOC_LAYOUT_PAULI); @@ -865,7 +839,6 @@ TEST_F(DFTUBaseTest, CompactWriterReaderCollinearRoundtrip) DFTU_BASE::write_occup_m(source, ucell, output, - 2, DFTU_BASE::SOC_LAYOUT_PAULI); } @@ -915,7 +888,6 @@ TEST_F(DFTUBaseTest, CompactWriterReaderSocPauliRoundtrip) DFTU_BASE::write_occup_m(source, ucell, output, - 4, DFTU_BASE::SOC_LAYOUT_PAULI); } @@ -964,7 +936,6 @@ TEST_F(DFTUBaseTest, CompactWriterReaderSocRealSpinBasisRoundtrip) DFTU_BASE::write_occup_m(source, ucell, output, - 4, DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); } @@ -1458,6 +1429,7 @@ TEST_F(DFTUBaseTest, AppendSnapshotGuardDftuOffBeforeIterValidation) // out_occ_mat=true but dft_plus_u=0: the guard must catch dft_plus_u // and return before the iter < 1 check fires. const DFTU_BASE::OccmatOutputCfg cfg = {1, 1, 5, true, 0}; + const DFTU_BASE::OccmatStepInfo step = {0, 0, 1e-6, 0.5}; // Use a unique subdirectory so files from other tests in the same // binary do not cause false positives. @@ -1468,14 +1440,9 @@ TEST_F(DFTUBaseTest, AppendSnapshotGuardDftuOffBeforeIterValidation) DFTU_BASE::append_ion_step_snapshot(dftu, ucell, out_dir, - 2, // nspin - 1, // npol - 0, // istep - 0, // iter — invalid, but DFT+U is off + step, false, false, - 1e-6, - 0.5, cfg, DFTU_BASE::SOC_LAYOUT_PAULI); From f0f02076fa1310b65261b661c7fbb0b40cfca841 Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Wed, 21 Oct 2026 23:20:25 +0800 Subject: [PATCH 5/6] Resolve CMake conflicts with develop formatting --- .../module_dftu/unittests/CMakeLists.txt | 347 +++++++++++------- source/source_pw/module_pwdft/CMakeLists.txt | 177 +++++---- .../module_pwdft/test/CMakeLists.txt | 204 +++++----- 3 files changed, 411 insertions(+), 317 deletions(-) diff --git a/source/source_lcao/module_dftu/unittests/CMakeLists.txt b/source/source_lcao/module_dftu/unittests/CMakeLists.txt index fb94f80174..437043ef60 100644 --- a/source/source_lcao/module_dftu/unittests/CMakeLists.txt +++ b/source/source_lcao/module_dftu/unittests/CMakeLists.txt @@ -7,154 +7,229 @@ abacus_disable_feature_definitions(__ROCM) abacus_disable_feature_definitions(__MLALGO) abacus_disable_feature_definitions(_OPENMP) -if (ENABLE_MPI) +if(ENABLE_MPI) -AddTest( - TARGET MODULE_DFTU_pots - LIBS base device - SOURCES test_dftu_nao_pots.cpp -) + addtest( + TARGET + MODULE_DFTU_pots + LIBS + base + device + SOURCES + test_dftu_nao_pots.cpp) -AddTest( - TARGET MODULE_DFTU_energy - LIBS base device - SOURCES test_dftu_nao_energy.cpp -) + addtest( + TARGET + MODULE_DFTU_energy + LIBS + base + device + SOURCES + test_dftu_nao_energy.cpp) -AddTest( - TARGET MODULE_DFTU_op_legacy - LIBS base device - SOURCES test_dftu_nao_op_legacy.cpp -) + addtest( + TARGET + MODULE_DFTU_op_legacy + LIBS + base + device + SOURCES + test_dftu_nao_op_legacy.cpp) -AddTest( - TARGET MODULE_DFTU_for_r - LIBS base device - SOURCES test_dftu_nao_for_r.cpp -) + addtest( + TARGET + MODULE_DFTU_for_r + LIBS + base + device + SOURCES + test_dftu_nao_for_r.cpp) -AddTest( - TARGET MODULE_DFTU_str_r - LIBS base device - SOURCES test_dftu_nao_str_r.cpp -) + addtest( + TARGET + MODULE_DFTU_str_r + LIBS + base + device + SOURCES + test_dftu_nao_str_r.cpp) -AddTest( - TARGET MODULE_DFTU_folding - LIBS base device - SOURCES test_dftu_nao_folding.cpp ../../../source_basis/module_ao/parallel_orbitals.cpp - # The test instantiates Parallel_Orbitals (a Parallel_2D subclass) and links - # the MPI-built base library, so it must keep __MPI to match the layout of - # Parallel_2D / Parallel_Orbitals across TUs (otherwise ODR/ABI mismatch). - KEEP_FEATURE_DEFINITIONS __MPI -) + addtest( + TARGET + MODULE_DFTU_folding + LIBS + base + device + SOURCES + test_dftu_nao_folding.cpp + ../../../source_basis/module_ao/parallel_orbitals.cpp + # The test instantiates Parallel_Orbitals (a Parallel_2D subclass) and links + # the MPI-built base library, so it must keep __MPI to match the layout of + # Parallel_2D / Parallel_Orbitals across TUs (otherwise ODR/ABI mismatch). + KEEP_FEATURE_DEFINITIONS + __MPI) -AddTest( - TARGET MODULE_DFTU_fs_r - LIBS base device - SOURCES test_dftu_nao_fs_r.cpp ../dftu_nao_fs_reduce.cpp -) + addtest( + TARGET + MODULE_DFTU_fs_r + LIBS + base + device + SOURCES + test_dftu_nao_fs_r.cpp + ../dftu_nao_fs_reduce.cpp) -AddTest( - TARGET MODULE_DFTU_fs_reduce - LIBS base device - SOURCES test_dftu_nao_fs_reduce.cpp ../dftu_nao_fs_reduce.cpp -) + addtest( + TARGET + MODULE_DFTU_fs_reduce + LIBS + base + device + SOURCES + test_dftu_nao_fs_reduce.cpp + ../dftu_nao_fs_reduce.cpp) -AddTest( - TARGET MODULE_DFTU_fs_k - LIBS base device - SOURCES test_dftu_nao_fs_k.cpp -) + addtest( + TARGET + MODULE_DFTU_fs_k + LIBS + base + device + SOURCES + test_dftu_nao_fs_k.cpp) -AddTest( - TARGET MODULE_DFTU_fs_accum - LIBS base device - SOURCES test_dftu_nao_fs_accum.cpp ../../../source_basis/module_ao/parallel_orbitals.cpp - # This test uses Parallel_2D::set_serial which is a real implementation - # in parallel_2d.cpp; keep __MPI so the serial path links correctly. - KEEP_FEATURE_DEFINITIONS __MPI -) + addtest( + TARGET + MODULE_DFTU_fs_accum + LIBS + base + device + SOURCES + test_dftu_nao_fs_accum.cpp + ../../../source_basis/module_ao/parallel_orbitals.cpp + # This test uses Parallel_2D::set_serial which is a real implementation + # in parallel_2d.cpp; keep __MPI so the serial path links correctly. + KEEP_FEATURE_DEFINITIONS + __MPI) -if(ENABLE_LCAO) -AddTest( - TARGET MODULE_DFTU_op - LIBS parameter psi base device container symmetry - SOURCES test_dftu_nao_op.cpp ../dftu_nao_op.cpp ../dftu_nao_adj.cpp ../dftu_nao_pots.cpp ../dftu_nao_fs_r.cpp ../dftu_nao_fs_reduce.cpp ../dftu_nao_for_r.cpp ../dftu_nao_str_r.cpp ../../../source_estate/module_dm/density_matrix.cpp - ../../../source_estate/module_dm/dmr_gamma.cpp - ../../../source_estate/module_dm/dmr_init.cpp - ../../../source_estate/module_dm/dm_setter.cpp - ../../../source_estate/module_dm/dm_getter.cpp - ../../../source_estate/module_dm/dm_tools.cpp - ../../../source_estate/module_dm/dmr_k.cpp - ../../../source_estate/module_dm/dmr_td.cpp - ../../../source_estate/module_dm/dmr_full.cpp - ../../../source_pw/module_pwdft/dftu_base.cpp - ../../../source_pw/module_pwdft/dftu_base_io.cpp - ../../../source_pw/module_pwdft/dftu_base_io_read.cpp - ../../../source_pw/module_pwdft/dftu_base_io_write.cpp - ../../../source_pw/module_pwdft/yukawa_screening.cpp - ../../../source_estate/occ_matrix.cpp - ../../../source_estate/occ_mixer.cpp - ../../../source_hamilt/module_hcontainer/func_folding.cpp - ../../../source_hamilt/module_hcontainer/base_matrix.cpp - ../../../source_hamilt/module_hcontainer/hcontainer.cpp - ../../../source_hamilt/module_hcontainer/atom_pair.cpp - ../../../source_basis/module_ao/parallel_orbitals.cpp - ../../../source_basis/module_ao/orb_atomic_lm.cpp - ../../module_operator_lcao/test/tmp_mocks.cpp - ../../../source_hamilt/operator.cpp - # This test calls MPI_Init in main() and drives ScaLAPACK distribution, - # so it must keep __MPI even though this directory disables it. - KEEP_FEATURE_DEFINITIONS __MPI -) + if(ENABLE_LCAO) + addtest( + TARGET + MODULE_DFTU_op + LIBS + parameter + psi + base + device + container + symmetry + SOURCES + test_dftu_nao_op.cpp + ../dftu_nao_op.cpp + ../dftu_nao_adj.cpp + ../dftu_nao_pots.cpp + ../dftu_nao_fs_r.cpp + ../dftu_nao_fs_reduce.cpp + ../dftu_nao_for_r.cpp + ../dftu_nao_str_r.cpp + ../../../source_estate/module_dm/density_matrix.cpp + ../../../source_estate/module_dm/dmr_gamma.cpp + ../../../source_estate/module_dm/dmr_init.cpp + ../../../source_estate/module_dm/dm_setter.cpp + ../../../source_estate/module_dm/dm_getter.cpp + ../../../source_estate/module_dm/dm_tools.cpp + ../../../source_estate/module_dm/dmr_k.cpp + ../../../source_estate/module_dm/dmr_td.cpp + ../../../source_estate/module_dm/dmr_full.cpp + ../../../source_pw/module_pwdft/dftu_base.cpp + ../../../source_pw/module_pwdft/dftu_base_io.cpp + ../../../source_pw/module_pwdft/dftu_base_io_read.cpp + ../../../source_pw/module_pwdft/dftu_base_io_write.cpp + ../../../source_pw/module_pwdft/yukawa_screening.cpp + ../../../source_estate/occ_matrix.cpp + ../../../source_estate/occ_mixer.cpp + ../../../source_hamilt/module_hcontainer/func_folding.cpp + ../../../source_hamilt/module_hcontainer/base_matrix.cpp + ../../../source_hamilt/module_hcontainer/hcontainer.cpp + ../../../source_hamilt/module_hcontainer/atom_pair.cpp + ../../../source_basis/module_ao/parallel_orbitals.cpp + ../../../source_basis/module_ao/orb_atomic_lm.cpp + ../../module_operator_lcao/test/tmp_mocks.cpp + ../../../source_hamilt/operator.cpp + # This test calls MPI_Init in main() and drives ScaLAPACK distribution, + # so it must keep __MPI even though this directory disables it. + KEEP_FEATURE_DEFINITIONS + __MPI) -AddTest( - TARGET MODULE_DFTU_ijr - LIBS parameter psi base device container - SOURCES test_dftu_nao_ijr.cpp - ../../../source_hamilt/module_hcontainer/func_folding.cpp - ../../../source_hamilt/module_hcontainer/base_matrix.cpp - ../../../source_hamilt/module_hcontainer/hcontainer.cpp - ../../../source_hamilt/module_hcontainer/atom_pair.cpp - ../../../source_basis/module_ao/parallel_orbitals.cpp - ../../../source_basis/module_ao/orb_atomic_lm.cpp - ../../module_operator_lcao/test/tmp_mocks.cpp - ../../../source_hamilt/operator.cpp - KEEP_FEATURE_DEFINITIONS __MPI -) + addtest( + TARGET + MODULE_DFTU_ijr + LIBS + parameter + psi + base + device + container + SOURCES + test_dftu_nao_ijr.cpp + ../../../source_hamilt/module_hcontainer/func_folding.cpp + ../../../source_hamilt/module_hcontainer/base_matrix.cpp + ../../../source_hamilt/module_hcontainer/hcontainer.cpp + ../../../source_hamilt/module_hcontainer/atom_pair.cpp + ../../../source_basis/module_ao/parallel_orbitals.cpp + ../../../source_basis/module_ao/orb_atomic_lm.cpp + ../../module_operator_lcao/test/tmp_mocks.cpp + ../../../source_hamilt/operator.cpp + KEEP_FEATURE_DEFINITIONS + __MPI) -AddTest( - TARGET MODULE_DFTU_base_io - LIBS parameter psi base device container symmetry - SOURCES test_dftu_base_io.cpp ../dftu_nao_op.cpp ../dftu_nao_adj.cpp ../dftu_nao_pots.cpp ../dftu_nao_fs_r.cpp ../dftu_nao_fs_reduce.cpp ../dftu_nao_for_r.cpp ../dftu_nao_str_r.cpp ../../../source_estate/module_dm/density_matrix.cpp - ../../../source_estate/module_dm/dmr_gamma.cpp - ../../../source_estate/module_dm/dmr_init.cpp - ../../../source_estate/module_dm/dm_setter.cpp - ../../../source_estate/module_dm/dm_getter.cpp - ../../../source_estate/module_dm/dm_tools.cpp - ../../../source_estate/module_dm/dmr_k.cpp - ../../../source_estate/module_dm/dmr_td.cpp - ../../../source_estate/module_dm/dmr_full.cpp - ../../../source_pw/module_pwdft/dftu_base.cpp - ../../../source_pw/module_pwdft/dftu_base_io.cpp - ../../../source_pw/module_pwdft/dftu_base_io_read.cpp - ../../../source_pw/module_pwdft/dftu_base_io_write.cpp - ../../../source_pw/module_pwdft/yukawa_screening.cpp - ../../../source_estate/occ_matrix.cpp - ../../../source_estate/occ_mixer.cpp - ../../../source_hamilt/module_hcontainer/func_folding.cpp - ../../../source_hamilt/module_hcontainer/base_matrix.cpp - ../../../source_hamilt/module_hcontainer/hcontainer.cpp - ../../../source_hamilt/module_hcontainer/atom_pair.cpp - ../../../source_basis/module_ao/parallel_orbitals.cpp - ../../../source_basis/module_ao/orb_atomic_lm.cpp - ../../module_operator_lcao/test/tmp_mocks.cpp - ../../../source_hamilt/operator.cpp - # Reuse the proven MODULE_DFTU_op link closure; MPI-built base/parallel - # objects must see __MPI across TUs. - KEEP_FEATURE_DEFINITIONS __MPI -) -endif() + addtest( + TARGET + MODULE_DFTU_base_io + LIBS + parameter + psi + base + device + container + symmetry + SOURCES + test_dftu_base_io.cpp + ../dftu_nao_op.cpp + ../dftu_nao_adj.cpp + ../dftu_nao_pots.cpp + ../dftu_nao_fs_r.cpp + ../dftu_nao_fs_reduce.cpp + ../dftu_nao_for_r.cpp + ../dftu_nao_str_r.cpp + ../../../source_estate/module_dm/density_matrix.cpp + ../../../source_estate/module_dm/dmr_gamma.cpp + ../../../source_estate/module_dm/dmr_init.cpp + ../../../source_estate/module_dm/dm_setter.cpp + ../../../source_estate/module_dm/dm_getter.cpp + ../../../source_estate/module_dm/dm_tools.cpp + ../../../source_estate/module_dm/dmr_k.cpp + ../../../source_estate/module_dm/dmr_td.cpp + ../../../source_estate/module_dm/dmr_full.cpp + ../../../source_pw/module_pwdft/dftu_base.cpp + ../../../source_pw/module_pwdft/dftu_base_io.cpp + ../../../source_pw/module_pwdft/dftu_base_io_read.cpp + ../../../source_pw/module_pwdft/dftu_base_io_write.cpp + ../../../source_pw/module_pwdft/yukawa_screening.cpp + ../../../source_estate/occ_matrix.cpp + ../../../source_estate/occ_mixer.cpp + ../../../source_hamilt/module_hcontainer/func_folding.cpp + ../../../source_hamilt/module_hcontainer/base_matrix.cpp + ../../../source_hamilt/module_hcontainer/hcontainer.cpp + ../../../source_hamilt/module_hcontainer/atom_pair.cpp + ../../../source_basis/module_ao/parallel_orbitals.cpp + ../../../source_basis/module_ao/orb_atomic_lm.cpp + ../../module_operator_lcao/test/tmp_mocks.cpp + ../../../source_hamilt/operator.cpp + # Reuse the proven MODULE_DFTU_op link closure; MPI-built base/parallel + # objects must see __MPI across TUs. + KEEP_FEATURE_DEFINITIONS + __MPI) + endif() endif() diff --git a/source/source_pw/module_pwdft/CMakeLists.txt b/source/source_pw/module_pwdft/CMakeLists.txt index b78a019700..2b0af1ec74 100644 --- a/source/source_pw/module_pwdft/CMakeLists.txt +++ b/source/source_pw/module_pwdft/CMakeLists.txt @@ -1,104 +1,101 @@ - -list(APPEND objects - hamilt_pw.cpp - pw_hybrid_nscf.cpp - pw_hybrid_source.cpp - td_history_pw.cpp - td_pw.cpp - op_pw_ekin.cpp - op_pw_ekin_td.cpp - op_pw_veff.cpp - projector_gradient.cpp - velocity_workspace.cpp - nonlocal_workspace.cpp - op_pw_nl.cpp - op_pw_nl_td.cpp - op_pw_meta.cpp - op_pw_vel.cpp - op_pw.cpp - op_pw_proj.cpp - op_pw_exx.cpp - op_pw_exx_ace.cpp - op_pw_exx_pot.cpp - kernels/exx_batch_op.cpp -) +list( + APPEND + objects + hamilt_pw.cpp + pw_hybrid_nscf.cpp + pw_hybrid_source.cpp + td_history_pw.cpp + td_pw.cpp + op_pw_ekin.cpp + op_pw_ekin_td.cpp + op_pw_veff.cpp + projector_gradient.cpp + velocity_workspace.cpp + nonlocal_workspace.cpp + op_pw_nl.cpp + op_pw_nl_td.cpp + op_pw_meta.cpp + op_pw_vel.cpp + op_pw.cpp + op_pw_proj.cpp + op_pw_exx.cpp + op_pw_exx_ace.cpp + op_pw_exx_pot.cpp + kernels/exx_batch_op.cpp) if(ENABLE_FLOAT_FFTW) - list(APPEND objects kernels/exx_batch_op_float.cpp) + list(APPEND objects kernels/exx_batch_op_float.cpp) else() - list(APPEND objects kernels/exx_batch_op_float_stub.cpp) + list(APPEND objects kernels/exx_batch_op_float_stub.cpp) endif() -list(APPEND objects - dftu_base.cpp - dftu_base_io.cpp - dftu_base_io_read.cpp - dftu_base_io_write.cpp - dftu_pw.cpp - dftu_pw_tools.cpp - yukawa_screening.cpp - setup_pot.cpp - setup_pwrho.cpp - setup_pwwfc.cpp - uspp_support.cpp - update_cell_pw.cpp - setup_dftu_pw.cpp - deltaspin_pw.cpp - force_pw_nl.cpp - force_pw_cc.cpp - force_pw_scc.cpp - force_pw.cpp - force_pw_us.cpp - force_pw_onsite.cpp - stress_cc.cpp - stress_ewa.cpp - stress_gga.cpp - stress_mgga.cpp - stress_har.cpp - stress_kin.cpp - stress_loc.cpp - stress_nl.cpp - stress_us.cpp - stress_onsite.cpp - stress_exx.cpp - stress_pw.cpp - vl_pw.cpp - vnl_pw.cpp - vnl_pw_alpha.cpp - vnl_pw_deeq.cpp - vnl_pw_getvnl.cpp - vnl_pw_grad.cpp - vnl_pw_init_vnl.cpp - vnl_pw_qrad.cpp +list( + APPEND + objects + dftu_base.cpp + dftu_base_io.cpp + dftu_base_io_read.cpp + dftu_base_io_write.cpp + dftu_pw.cpp + dftu_pw_tools.cpp + yukawa_screening.cpp + setup_pot.cpp + setup_pwrho.cpp + setup_pwwfc.cpp + uspp_support.cpp + update_cell_pw.cpp + setup_dftu_pw.cpp + deltaspin_pw.cpp + force_pw_nl.cpp + force_pw_cc.cpp + force_pw_scc.cpp + force_pw.cpp + force_pw_us.cpp + force_pw_onsite.cpp + stress_cc.cpp + stress_ewa.cpp + stress_gga.cpp + stress_mgga.cpp + stress_har.cpp + stress_kin.cpp + stress_loc.cpp + stress_nl.cpp + stress_us.cpp + stress_onsite.cpp + stress_exx.cpp + stress_pw.cpp + vl_pw.cpp + vnl_pw.cpp + vnl_pw_alpha.cpp + vnl_pw_deeq.cpp + vnl_pw_getvnl.cpp + vnl_pw_grad.cpp + vnl_pw_init_vnl.cpp + vnl_pw_qrad.cpp vnl_pw_table.cpp - stru_fac.cpp - stru_fac_k.cpp - soc.cpp - elecond.cpp - fs_nonlocal_tools.cpp - fs_nonlocal_tools_vkb.cpp - fs_nonlocal_tools_stress.cpp - fs_nonlocal_tools_force.cpp - fs_kin_tools.cpp - nonlocal_maths.cpp - onsite_proj_force_stress.cpp - vsep_pw.cpp -) + stru_fac.cpp + stru_fac_k.cpp + soc.cpp + elecond.cpp + fs_nonlocal_tools.cpp + fs_nonlocal_tools_vkb.cpp + fs_nonlocal_tools_stress.cpp + fs_nonlocal_tools_force.cpp + fs_kin_tools.cpp + nonlocal_maths.cpp + onsite_proj_force_stress.cpp + vsep_pw.cpp) -add_library( - module_pwdft - OBJECT - ${objects} -) +add_library(module_pwdft OBJECT ${objects}) if(ENABLE_COVERAGE) add_coverage(module_pwdft) endif() -if (BUILD_TESTING) - if(ENABLE_MPI) - add_subdirectory(test) - add_subdirectory(kernels/test) - endif() +if(BUILD_TESTING) + if(ENABLE_MPI) + add_subdirectory(test) + add_subdirectory(kernels/test) + endif() endif() diff --git a/source/source_pw/module_pwdft/test/CMakeLists.txt b/source/source_pw/module_pwdft/test/CMakeLists.txt index 1b0c556dc7..91e0ce07bc 100644 --- a/source/source_pw/module_pwdft/test/CMakeLists.txt +++ b/source/source_pw/module_pwdft/test/CMakeLists.txt @@ -3,100 +3,122 @@ abacus_disable_feature_definitions(__CUDA) abacus_disable_feature_definitions(__ROCM) abacus_disable_feature_definitions(__EXX) -AddTest( - TARGET MODULE_PW_pwdft_soc - LIBS parameter - SOURCES soc_test.cpp ../soc.cpp - ../../../source_base/global_variable.cpp - ../../../source_base/global_function.cpp - ../../../source_base/global_file.cpp - ../../../source_base/tool_quit.cpp - ../../../source_base/global_variable.cpp - ../../../source_base/global_file.cpp - ../../../source_base/memory_recorder.cpp - ../../../source_base/timer.cpp - ../../../source_base/module_external/blas_connector_base.cpp ../../../source_base/module_external/blas_connector_vector.cpp ../../../source_base/module_external/blas_connector_matrix.cpp - ../../../source_base/parallel_global.cpp - ../../../source_base/parallel_comm.cpp - ../../../source_base/parallel_common.cpp - ../../../source_base/parallel_reduce.cpp -) +addtest( + TARGET + MODULE_PW_pwdft_soc + LIBS + parameter + SOURCES + soc_test.cpp + ../soc.cpp + ../../../source_base/global_variable.cpp + ../../../source_base/global_function.cpp + ../../../source_base/global_file.cpp + ../../../source_base/tool_quit.cpp + ../../../source_base/global_variable.cpp + ../../../source_base/global_file.cpp + ../../../source_base/memory_recorder.cpp + ../../../source_base/timer.cpp + ../../../source_base/module_external/blas_connector_base.cpp + ../../../source_base/module_external/blas_connector_vector.cpp + ../../../source_base/module_external/blas_connector_matrix.cpp + ../../../source_base/parallel_global.cpp + ../../../source_base/parallel_comm.cpp + ../../../source_base/parallel_common.cpp + ../../../source_base/parallel_reduce.cpp) -AddTest( - TARGET MODULE_PW_dftu_base_test - LIBS parameter base device - SOURCES dftu_base_test.cpp - ../dftu_base.cpp - ../dftu_base_io.cpp - ../dftu_base_io_read.cpp - ../dftu_base_io_write.cpp - ../dftu_pw_tools.cpp - ../yukawa_screening.cpp - ../../../source_estate/occ_matrix.cpp - ../../../source_estate/occ_mixer.cpp - ../../../source_cell/unitcell.cpp - ../../../source_cell/unitcell_index.cpp - ../../../source_cell/unitcell_stats.cpp - ../../../source_cell/unitcell_setup.cpp - ../../../source_cell/atom_spec.cpp - ../../../source_cell/atom_pseudo.cpp - ../../../source_cell/pseudo.cpp - ../../../source_cell/magnetism.cpp - ../../../source_cell/sep.cpp - ../../../source_cell/sep_cell.cpp - ../../../source_cell/read_atom_species.cpp - ../../../source_cell/read_atoms.cpp - ../../../source_cell/read_atoms_helper.cpp - ../../../source_cell/cell_tools.cpp - ../../../source_cell/print_cell.cpp - ../../../source_cell/read_orb.cpp - ../../../source_cell/read_stru.cpp - ../../../source_cell/bcast_cell.cpp -) - -AddTest( - TARGET MODULE_PW_test_dftu_pw_tools - LIBS parameter base device - SOURCES test_dftu_pw_tools.cpp - ../dftu_pw_tools.cpp -) +addtest( + TARGET + MODULE_PW_dftu_base_test + LIBS + parameter + base + device + SOURCES + dftu_base_test.cpp + ../dftu_base.cpp + ../dftu_base_io.cpp + ../dftu_base_io_read.cpp + ../dftu_base_io_write.cpp + ../dftu_pw_tools.cpp + ../yukawa_screening.cpp + ../../../source_estate/occ_matrix.cpp + ../../../source_estate/occ_mixer.cpp + ../../../source_cell/unitcell.cpp + ../../../source_cell/unitcell_index.cpp + ../../../source_cell/unitcell_stats.cpp + ../../../source_cell/unitcell_setup.cpp + ../../../source_cell/atom_spec.cpp + ../../../source_cell/atom_pseudo.cpp + ../../../source_cell/pseudo.cpp + ../../../source_cell/magnetism.cpp + ../../../source_cell/sep.cpp + ../../../source_cell/sep_cell.cpp + ../../../source_cell/read_atom_species.cpp + ../../../source_cell/read_atoms.cpp + ../../../source_cell/read_atoms_helper.cpp + ../../../source_cell/cell_tools.cpp + ../../../source_cell/print_cell.cpp + ../../../source_cell/read_orb.cpp + ../../../source_cell/read_stru.cpp + ../../../source_cell/bcast_cell.cpp) -AddTest( - TARGET MODULE_PW_structure_factor_test - LIBS parameter base device planewave - SOURCES structure_factor_test.cpp ../stru_fac.cpp - ../../../source_cell/unitcell.cpp - ../../../source_cell/unitcell_index.cpp - ../../../source_cell/unitcell_stats.cpp - ../../../source_cell/unitcell_setup.cpp +addtest( + TARGET + MODULE_PW_test_dftu_pw_tools + LIBS + parameter + base + device + SOURCES + test_dftu_pw_tools.cpp + ../dftu_pw_tools.cpp) - ../../../source_cell/update_cell.cpp - ../../../source_cell/bcast_cell.cpp - ../../../source_cell/print_cell.cpp - ../../../source_cell/atom_spec.cpp - ../../../source_cell/atom_pseudo.cpp - ../../../source_cell/pseudo.cpp - ../../../source_cell/read_stru.cpp - ../../../source_cell/read_atom_species.cpp - ../../../source_cell/read_atoms.cpp - ../../../source_cell/read_atoms_helper.cpp - ../../../source_cell/cell_tools.cpp - ../../../source_cell/read_pp.cpp - ../../../source_cell/read_pp_complete.cpp - ../../../source_cell/read_pp_upf100.cpp - ../../../source_cell/read_pp_upf201.cpp - ../../../source_cell/read_pp_vwr.cpp - ../../../source_cell/read_pp_blps.cpp +addtest( + TARGET + MODULE_PW_structure_factor_test + LIBS + parameter + base + device + planewave + SOURCES + structure_factor_test.cpp + ../stru_fac.cpp + ../../../source_cell/unitcell.cpp + ../../../source_cell/unitcell_index.cpp + ../../../source_cell/unitcell_stats.cpp + ../../../source_cell/unitcell_setup.cpp + ../../../source_cell/update_cell.cpp + ../../../source_cell/bcast_cell.cpp + ../../../source_cell/print_cell.cpp + ../../../source_cell/atom_spec.cpp + ../../../source_cell/atom_pseudo.cpp + ../../../source_cell/pseudo.cpp + ../../../source_cell/read_stru.cpp + ../../../source_cell/read_atom_species.cpp + ../../../source_cell/read_atoms.cpp + ../../../source_cell/read_atoms_helper.cpp + ../../../source_cell/cell_tools.cpp + ../../../source_cell/read_pp.cpp + ../../../source_cell/read_pp_complete.cpp + ../../../source_cell/read_pp_upf100.cpp + ../../../source_cell/read_pp_upf201.cpp + ../../../source_cell/read_pp_vwr.cpp + ../../../source_cell/read_pp_blps.cpp ../../../source_cell/sep.cpp ../../../source_cell/sep_cell.cpp - ../../../source_cell/read_pp_ucell.cpp - ../../../source_cell/cal_wfc.cpp - ../../../source_cell/cal_nelec_nband.cpp - ../../../source_cell/read_orb.cpp -) + ../../../source_cell/read_pp_ucell.cpp + ../../../source_cell/cal_wfc.cpp + ../../../source_cell/cal_nelec_nband.cpp + ../../../source_cell/read_orb.cpp) -AddTest( - TARGET MODULE_PW_pw_hybrid_source - LIBS base device - SOURCES test_pw_hybrid_source.cpp ../pw_hybrid_source.cpp -) +addtest( + TARGET + MODULE_PW_pw_hybrid_source + LIBS + base + device + SOURCES + test_pw_hybrid_source.cpp + ../pw_hybrid_source.cpp) From 98cd90018d026030b5e8f6c81b8361d3ac54dbea Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Thu, 22 Oct 2026 07:14:23 +0800 Subject: [PATCH 6/6] Address review: assert mag/nelec diagnostics, drop dead code - Remove the trailing first_scf_step return and its output-cadence comment from iter_init_dftu_pw; output() was removed earlier, so the guard no longer controls anything. - Extend the compact writer roundtrip tests to assert the serialized mag and nelec header values against values computed from the source occupation matrix, including the nonzero y magnetization in the SOC Pauli case. read_occup_m only checks these fields for finiteness, so matrix-only comparisons could not detect an incorrect diagnostic. --- .../source_pw/module_pwdft/setup_dftu_pw.cpp | 7 -- .../module_pwdft/test/dftu_base_test.cpp | 88 +++++++++++++++++++ 2 files changed, 88 insertions(+), 7 deletions(-) diff --git a/source/source_pw/module_pwdft/setup_dftu_pw.cpp b/source/source_pw/module_pwdft/setup_dftu_pw.cpp index d3be1285fa..877445b88b 100644 --- a/source/source_pw/module_pwdft/setup_dftu_pw.cpp +++ b/source/source_pw/module_pwdft/setup_dftu_pw.cpp @@ -66,13 +66,6 @@ void iter_init_dftu_pw(const int iter, dftu.get_uterm_mat(), dftu.energy_ref()); } - // Keep the historical output cadence: output() is skipped at the very - // first electronic step. - if (first_scf_step) - { - return; - } - } } diff --git a/source/source_pw/module_pwdft/test/dftu_base_test.cpp b/source/source_pw/module_pwdft/test/dftu_base_test.cpp index decf278af8..596da0b7ca 100644 --- a/source/source_pw/module_pwdft/test/dftu_base_test.cpp +++ b/source/source_pw/module_pwdft/test/dftu_base_test.cpp @@ -865,6 +865,46 @@ TEST_F(DFTUBaseTest, CompactWriterReaderCollinearRoundtrip) } } } + + // Assert the diagnostic header values that read_occup_m only validates + // for finiteness: per-spin nelec (trace) and collinear mag (up - down). + double expected_nelec[2] = {0.0, 0.0}; + for (int spin = 0; spin < 2; ++spin) + { + for (int m = 0; m < 5; ++m) + { + expected_nelec[spin] += source.occmat().get(0, 2, spin, m, m); + } + } + const double expected_mag = expected_nelec[0] - expected_nelec[1]; + + std::ifstream input(filename.c_str()); + ASSERT_TRUE(input.is_open()); + std::string token; + double parsed_mag = 0.0; + double parsed_nelec[2] = {0.0, 0.0}; + int spin_index = -1; + while (input >> token) + { + if (token == "mag") + { + ASSERT_TRUE(input >> parsed_mag); + } + else if (token == "spin") + { + ASSERT_TRUE(input >> spin_index); + ASSERT_GE(spin_index, 1); + ASSERT_LE(spin_index, 2); + } + else if (token == "nelec") + { + ASSERT_EQ(spin_index == 1 || spin_index == 2, true); + ASSERT_TRUE(input >> parsed_nelec[spin_index - 1]); + } + } + EXPECT_NEAR(parsed_mag, expected_mag, 1e-7); + EXPECT_NEAR(parsed_nelec[0], expected_nelec[0], 1e-7); + EXPECT_NEAR(parsed_nelec[1], expected_nelec[1], 1e-7); } TEST_F(DFTUBaseTest, CompactWriterReaderSocPauliRoundtrip) @@ -908,6 +948,54 @@ TEST_F(DFTUBaseTest, CompactWriterReaderSocPauliRoundtrip) source.occmat().mat(0, 2, 0).c[index], 1e-8); } + + // Assert the SOC diagnostic header: mag x/y/z and per-spin nelec. + // Pauli blocks are charge/spin_x/spin_y/spin_z, so the writer reports + // mag_x = sum(spin_x diag), mag_y = sum(spin_y diag), mag_z = sum(spin_z diag) + // and nelec = 0.5 * sum(charge +/- spin_z) on the diagonal. + const ModuleBase::matrix& occ = source.occmat().mat(0, 2, 0); + double expected_mag[3] = {0.0, 0.0, 0.0}; + double expected_nelec_soc[2] = {0.0, 0.0}; + for (int m = 0; m < 5; ++m) + { + const int diag = m * 5 + m; + expected_mag[0] += occ.c[1 * block_size + diag]; + expected_mag[1] += occ.c[2 * block_size + diag]; + expected_mag[2] += occ.c[3 * block_size + diag]; + expected_nelec_soc[0] += 0.5 * (occ.c[0 * block_size + diag] + occ.c[3 * block_size + diag]); + expected_nelec_soc[1] += 0.5 * (occ.c[0 * block_size + diag] - occ.c[3 * block_size + diag]); + } + // The fixture fills every block entry, including a nonzero spin_y block, + // so mag_y must be nonzero for this assertion to be meaningful. + EXPECT_NE(expected_mag[1], 0.0); + + std::ifstream soc_input(filename.c_str()); + ASSERT_TRUE(soc_input.is_open()); + std::string token; + double parsed_mag_soc[3] = {0.0, 0.0, 0.0}; + double parsed_nelec_soc[2] = {0.0, 0.0}; + int spin_index = -1; + while (soc_input >> token) + { + if (token == "mag") + { + ASSERT_TRUE(soc_input >> parsed_mag_soc[0] >> parsed_mag_soc[1] >> parsed_mag_soc[2]); + } + else if (token == "spin") + { + ASSERT_TRUE(soc_input >> spin_index); + } + else if (token == "nelec") + { + ASSERT_EQ(spin_index == 1 || spin_index == 2, true); + ASSERT_TRUE(soc_input >> parsed_nelec_soc[spin_index - 1]); + } + } + EXPECT_NEAR(parsed_mag_soc[0], expected_mag[0], 1e-7); + EXPECT_NEAR(parsed_mag_soc[1], expected_mag[1], 1e-7); + EXPECT_NEAR(parsed_mag_soc[2], expected_mag[2], 1e-7); + EXPECT_NEAR(parsed_nelec_soc[0], expected_nelec_soc[0], 1e-7); + EXPECT_NEAR(parsed_nelec_soc[1], expected_nelec_soc[1], 1e-7); } TEST_F(DFTUBaseTest, CompactWriterReaderSocRealSpinBasisRoundtrip)