From b08549ea7b668705e48f9f30bf2891787590e712 Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Sat, 17 Oct 2026 20:30:22 +0800 Subject: [PATCH 01/13] Refactor DFT+U occupation reader into separate source --- .../module_dftu/unittests/CMakeLists.txt | 2 + source/source_pw/module_pwdft/CMakeLists.txt | 1 + .../source_pw/module_pwdft/dftu_base_io.cpp | 347 ----------------- .../module_pwdft/dftu_base_io_read.cpp | 367 ++++++++++++++++++ .../module_pwdft/test/CMakeLists.txt | 1 + 5 files changed, 371 insertions(+), 347 deletions(-) create mode 100644 source/source_pw/module_pwdft/dftu_base_io_read.cpp diff --git a/source/source_lcao/module_dftu/unittests/CMakeLists.txt b/source/source_lcao/module_dftu/unittests/CMakeLists.txt index fb7d40cc6c..6a50a09ffe 100644 --- a/source/source_lcao/module_dftu/unittests/CMakeLists.txt +++ b/source/source_lcao/module_dftu/unittests/CMakeLists.txt @@ -91,6 +91,7 @@ AddTest( ../../../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/yukawa_screening.cpp ../../../source_estate/occ_matrix.cpp ../../../source_estate/occ_mixer.cpp @@ -136,6 +137,7 @@ AddTest( ../../../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/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 cf652068f3..404fb05b4f 100644 --- a/source/source_pw/module_pwdft/CMakeLists.txt +++ b/source/source_pw/module_pwdft/CMakeLists.txt @@ -31,6 +31,7 @@ endif() list(APPEND objects dftu_base.cpp dftu_base_io.cpp + dftu_base_io_read.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 b0cbb9e320..05a13700ac 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.cpp +++ b/source/source_pw/module_pwdft/dftu_base_io.cpp @@ -227,41 +227,6 @@ void write_provenance_header(std::ostream& os, const int istep) os << "# RELAX STEP " << istep + 1 << "\n"; } -// Consume the label of one SOC block in the compact layout, e.g. -// "spin 1 nelec 1.07922327" or "spin 12 re". The optional -// "nelec " token pair is skipped when present. -void skip_compact_soc_label(std::ifstream& ifs, char* word) -{ - ifs >> word; // "spin" - ifs >> word; // block index: "1", "2" or "12" - ifs >> word; // "nelec", "re" or "im" - if (strcmp(word, "nelec") == 0) - { - double nelec = 0.0; - ifs >> nelec; - } -} - -// Read n*n doubles of one compact matrix block row by row. When out is -// non-null the values are stored, otherwise they are read and discarded. -void read_matrix_block(std::ifstream& ifs, const int n, double* out) -{ - for (int i = 0; i < n * n; ++i) - { - double value = 0.0; - ifs >> value; - if (ifs.fail()) - { - ModuleBase::WARNING_QUIT("read_matrix_block", - "failed to read a matrix element from the occupation matrix file"); - } - if (out != nullptr) - { - out[i] = value; - } - } -} - } // namespace @@ -308,318 +273,6 @@ std::string find_first_existing_file(const std::string& dir, return std::string(); } -void read_occup_m(const UnitCell& ucell, - OccupationMatrix& occ, - const std::vector& l_channel, - const int init_occ_mat, - const std::string& fn, - const std::string& init_chg, - int nspin, - int npol, - OccmatSocLayout soc_layout) -{ - ModuleBase::TITLE("DFTU_BASE", "read_occup_m"); - - if (GlobalV::MY_RANK != 0) - { - return; - } - - std::ifstream ifdftu(fn.c_str(), std::ios::in); - - if (!ifdftu) - { - if (init_occ_mat > 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "Can not find the file dm_onsite_ini.txt. Please check your dm_onsite_ini.txt"); - } - else - { - if (init_chg == "file") - { - const std::string not_found_msg = "Can not find the file " + fn - + ". Please do scf calculation first"; - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", not_found_msg); - } - } - const std::string open_failed_msg = "Can not open " + fn + " file"; - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", open_failed_msg); - } - - GlobalV::ofs_running << " DFT+U: read occupation matrix from " << fn << std::endl; - std::cout << " DFT+U: read occupation matrix from " << fn << std::endl; - - ifdftu.clear(); - ifdftu.seekg(0); - - char word[20]; - - int T = 0; - int iat = 0; - int spin = 0; - int L = 0; - int zeta = 0; - - while (true) - { - if (!(ifdftu >> word)) - { - break; - } - - // Comment lines ("# ...") carry the provenance header and the - // electronic-step marker in the compact layout; skip the line. - if (word[0] == '#') - { - ifdftu.ignore(150, '\n'); - continue; - } - - // The legacy "Atom=" token and the compact "Atom" label both - // start one atom block; every other token (element prefix, log - // sentences) is ignored. - const bool legacy_atom = (strcmp(word, "Atom=") == 0); - const bool compact_atom = (strcmp(word, "Atom") == 0); - if (!legacy_atom && !compact_atom) - { - continue; - } - const bool legacy_layout = legacy_atom; - - ifdftu >> iat; - iat -= 1; - - if (legacy_layout) - { - ifdftu >> word; - - // Accept both the compact "L=2" and the legacy split "L= 2". - if (strncmp("L=", word, 2) != 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); - } - if (strlen(word) > 2) - { - L = atoi(word + 2); - } - else - { - ifdftu >> L; - } - ifdftu >> word; - - // The ORBITAL= token is optional in newly written files; when - // absent the radial channel defaults to 0. - if (strncmp("ORBITAL=", word, 8) == 0) - { - if (strlen(word) > 8) - { - zeta = atoi(word + 8); - } - else - { - ifdftu >> zeta; - } - } - else - { - zeta = 0; - } - ifdftu.ignore(150, '\n'); - - if (zeta != 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", - "only the first radial channel (ORBITAL=0) is supported"); - } - } - else - { - // Compact header, e.g. " Fe Atom 1 L 2 mag 0.03..." - ifdftu >> word; - if (strcmp(word, "L") != 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); - } - ifdftu >> L; - // The optional "mag ..." tail is redundant for a restart. - ifdftu.ignore(150, '\n'); - zeta = 0; - } - - T = ucell.iat2it[iat]; - const int NL = ucell.atoms[T].nwl + 1; - - for (int l = 0; l < NL; l++) - { - if (l != l_channel[T]) - { - continue; - } - - const int nm = 2 * L + 1; - - if (nspin == 1 || nspin == 2) - { - for (int is = 0; is < 2; is++) - { - if (legacy_layout) - { - ifdftu >> word; - if (strcmp("spin=", word) != 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); - } - ifdftu >> spin; - spin -= 1; - ifdftu.ignore(150, '\n'); - - double value = 0.0; - for (int m0 = 0; m0 < nm; m0++) - { - for (int m1 = 0; m1 < nm; m1++) - { - ifdftu >> value; - occ.set(iat, L, spin, m0, m1, value); - } - ifdftu.ignore(150, '\n'); - } - } - else - { - // Compact label: "spin 1 nelec 0.468..." - ifdftu >> word; - if (strcmp("spin", word) != 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); - } - int spin_idx = 0; - ifdftu >> spin_idx; - if (spin_idx != is + 1) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG SPIN INDEX IN Plus_U FILE"); - } - ifdftu >> word; - if (strcmp("nelec", word) != 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); - } - double nelec = 0.0; - ifdftu >> nelec; - - double value = 0.0; - for (int m0 = 0; m0 < nm; m0++) - { - for (int m1 = 0; m1 < nm; m1++) - { - ifdftu >> value; - occ.set(iat, L, is, m0, m1, value); - } - } - } - } - } - else if (nspin == 4) // SOC - { - if (legacy_layout) - { - double value = 0.0; - for (int m0 = 0; m0 < nm; m0++) - { - for (int ipol0 = 0; ipol0 < npol; ipol0++) - { - const int m0_all = m0 + nm * ipol0; - - for (int m1 = 0; m1 < nm; m1++) - { - for (int ipol1 = 0; ipol1 < npol; ipol1++) - { - int m1_all = m1 + nm * ipol1; - ifdftu >> value; - occ.set(iat, L, 0, m0_all, m1_all, value); - } - } - ifdftu.ignore(150, '\n'); - } - } - } - else - { - // Four labeled blocks, in the same order they are - // written: up/up, up/down Re, up/down Im, down/down. - std::vector uu(nm * nm); - std::vector re(nm * nm); - std::vector im(nm * nm); - std::vector dd(nm * nm); - - skip_compact_soc_label(ifdftu, word); - read_matrix_block(ifdftu, nm, uu.data()); - - skip_compact_soc_label(ifdftu, word); - read_matrix_block(ifdftu, nm, re.data()); - - // The imaginary block is needed by the PW (Pauli) - // path to reconstruct b2 = 2*Im(n_ud); the LCAO - // (real spin-basis) path discards it because its - // storage can not represent Im(n_ud). - skip_compact_soc_label(ifdftu, word); - read_matrix_block(ifdftu, nm, im.data()); - - skip_compact_soc_label(ifdftu, word); - read_matrix_block(ifdftu, nm, dd.data()); - - if (soc_layout == SOC_LAYOUT_PAULI) - { - // Reconstruct the 4 contiguous Pauli blocks - // [b0, b1, b2, b3] in the 2m x 2m flat buffer - // from the file's (n_uu, Re(n_ud), Im(n_ud), - // n_dd) representation. The writer prints - // n_uu = (b0 + b3)/2 - // Re(n_ud) = b1/2 - // Im(n_ud) = b2/2 - // n_dd = (b0 - b3)/2 - // so the inverse is - // b0 = n_uu + n_dd - // b1 = 2 * Re(n_ud) - // b2 = 2 * Im(n_ud) - // b3 = n_uu - n_dd. - ModuleBase::matrix& occ0 = occ.mat(iat, L, 0); - const int m2 = nm * nm; - for (int k = 0; k < m2; ++k) - { - occ0.c[0 * m2 + k] = uu[k] + dd[k]; - occ0.c[1 * m2 + k] = 2.0 * re[k]; - occ0.c[2 * m2 + k] = 2.0 * im[k]; - occ0.c[3 * m2 + k] = uu[k] - dd[k]; - } - } - else // SOC_LAYOUT_SPIN_BASIS_REAL - { - // Reconstruct the real 2m x 2m spin-basis matrix. - // Mirror the legacy writer: the Re block fills both - // (up, down) and (down, up) at the same (m0, m1). - // Im(n_ud) is dropped because the storage can not - // represent it. - for (int m0 = 0; m0 < nm; m0++) - { - for (int m1 = 0; m1 < nm; m1++) - { - const int k = m0 * nm + m1; - occ.set(iat, L, 0, m0, m1, uu[k]); - occ.set(iat, L, 0, m0, nm + m1, re[k]); - occ.set(iat, L, 0, nm + m0, m1, re[k]); - occ.set(iat, L, 0, nm + m0, nm + m1, dd[k]); - } - } - } - } - } - } - } - - return; -} - #ifdef __MPI /// Broadcast the local occupation number matrices from rank 0 to all ranks. /// diff --git a/source/source_pw/module_pwdft/dftu_base_io_read.cpp b/source/source_pw/module_pwdft/dftu_base_io_read.cpp new file mode 100644 index 0000000000..a94d2ef2d8 --- /dev/null +++ b/source/source_pw/module_pwdft/dftu_base_io_read.cpp @@ -0,0 +1,367 @@ +#include "source_pw/module_pwdft/dftu_base_io.h" + +#include "source_base/global_function.h" +#include "source_base/global_variable.h" +#include "source_cell/unitcell.h" + +#include +#include +#include +#include + +namespace +{ + +// Consume the label of one SOC block in the compact layout, e.g. +// "spin 1 nelec 1.07922327" or "spin 12 re". The optional +// "nelec " token pair is skipped when present. +void skip_compact_soc_label(std::ifstream& ifs, char* word) +{ + ifs >> word; // "spin" + ifs >> word; // block index: "1", "2" or "12" + ifs >> word; // "nelec", "re" or "im" + if (strcmp(word, "nelec") == 0) + { + double nelec = 0.0; + ifs >> nelec; + } +} + +// Read n*n doubles of one compact matrix block row by row. When out is +// non-null the values are stored, otherwise they are read and discarded. +void read_matrix_block(std::ifstream& ifs, const int n, double* out) +{ + for (int i = 0; i < n * n; ++i) + { + double value = 0.0; + ifs >> value; + if (ifs.fail()) + { + ModuleBase::WARNING_QUIT("read_matrix_block", + "failed to read a matrix element from the occupation matrix file"); + } + if (out != nullptr) + { + out[i] = value; + } + } +} + +} // namespace + +namespace DFTU_BASE +{ + +void read_occup_m(const UnitCell& ucell, + OccupationMatrix& occ, + const std::vector& l_channel, + const int init_occ_mat, + const std::string& fn, + const std::string& init_chg, + int nspin, + int npol, + OccmatSocLayout soc_layout) +{ + ModuleBase::TITLE("DFTU_BASE", "read_occup_m"); + + if (GlobalV::MY_RANK != 0) + { + return; + } + + std::ifstream ifdftu(fn.c_str(), std::ios::in); + + if (!ifdftu) + { + if (init_occ_mat > 0) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "Can not find the file dm_onsite_ini.txt. Please check your dm_onsite_ini.txt"); + } + else + { + if (init_chg == "file") + { + const std::string not_found_msg = "Can not find the file " + fn + + ". Please do scf calculation first"; + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", not_found_msg); + } + } + const std::string open_failed_msg = "Can not open " + fn + " file"; + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", open_failed_msg); + } + + GlobalV::ofs_running << " DFT+U: read occupation matrix from " << fn << std::endl; + std::cout << " DFT+U: read occupation matrix from " << fn << std::endl; + + ifdftu.clear(); + ifdftu.seekg(0); + + char word[20]; + + int T = 0; + int iat = 0; + int spin = 0; + int L = 0; + int zeta = 0; + + while (true) + { + if (!(ifdftu >> word)) + { + break; + } + + // Comment lines ("# ...") carry the provenance header and the + // electronic-step marker in the compact layout; skip the line. + if (word[0] == '#') + { + ifdftu.ignore(150, '\n'); + continue; + } + + // The legacy "Atom=" token and the compact "Atom" label both + // start one atom block; every other token (element prefix, log + // sentences) is ignored. + const bool legacy_atom = (strcmp(word, "Atom=") == 0); + const bool compact_atom = (strcmp(word, "Atom") == 0); + if (!legacy_atom && !compact_atom) + { + continue; + } + const bool legacy_layout = legacy_atom; + + ifdftu >> iat; + iat -= 1; + + if (legacy_layout) + { + ifdftu >> word; + + // Accept both the compact "L=2" and the legacy split "L= 2". + if (strncmp("L=", word, 2) != 0) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); + } + if (strlen(word) > 2) + { + L = atoi(word + 2); + } + else + { + ifdftu >> L; + } + ifdftu >> word; + + // The ORBITAL= token is optional in newly written files; when + // absent the radial channel defaults to 0. + if (strncmp("ORBITAL=", word, 8) == 0) + { + if (strlen(word) > 8) + { + zeta = atoi(word + 8); + } + else + { + ifdftu >> zeta; + } + } + else + { + zeta = 0; + } + ifdftu.ignore(150, '\n'); + + if (zeta != 0) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", + "only the first radial channel (ORBITAL=0) is supported"); + } + } + else + { + // Compact header, e.g. " Fe Atom 1 L 2 mag 0.03..." + ifdftu >> word; + if (strcmp(word, "L") != 0) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); + } + ifdftu >> L; + // The optional "mag ..." tail is redundant for a restart. + ifdftu.ignore(150, '\n'); + zeta = 0; + } + + T = ucell.iat2it[iat]; + const int NL = ucell.atoms[T].nwl + 1; + + for (int l = 0; l < NL; l++) + { + if (l != l_channel[T]) + { + continue; + } + + const int nm = 2 * L + 1; + + if (nspin == 1 || nspin == 2) + { + for (int is = 0; is < 2; is++) + { + if (legacy_layout) + { + ifdftu >> word; + if (strcmp("spin=", word) != 0) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); + } + ifdftu >> spin; + spin -= 1; + ifdftu.ignore(150, '\n'); + + double value = 0.0; + for (int m0 = 0; m0 < nm; m0++) + { + for (int m1 = 0; m1 < nm; m1++) + { + ifdftu >> value; + occ.set(iat, L, spin, m0, m1, value); + } + ifdftu.ignore(150, '\n'); + } + } + else + { + // Compact label: "spin 1 nelec 0.468..." + ifdftu >> word; + if (strcmp("spin", word) != 0) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); + } + int spin_idx = 0; + ifdftu >> spin_idx; + if (spin_idx != is + 1) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG SPIN INDEX IN Plus_U FILE"); + } + ifdftu >> word; + if (strcmp("nelec", word) != 0) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); + } + double nelec = 0.0; + ifdftu >> nelec; + + double value = 0.0; + for (int m0 = 0; m0 < nm; m0++) + { + for (int m1 = 0; m1 < nm; m1++) + { + ifdftu >> value; + occ.set(iat, L, is, m0, m1, value); + } + } + } + } + } + else if (nspin == 4) // SOC + { + if (legacy_layout) + { + double value = 0.0; + for (int m0 = 0; m0 < nm; m0++) + { + for (int ipol0 = 0; ipol0 < npol; ipol0++) + { + const int m0_all = m0 + nm * ipol0; + + for (int m1 = 0; m1 < nm; m1++) + { + for (int ipol1 = 0; ipol1 < npol; ipol1++) + { + int m1_all = m1 + nm * ipol1; + ifdftu >> value; + occ.set(iat, L, 0, m0_all, m1_all, value); + } + } + ifdftu.ignore(150, '\n'); + } + } + } + else + { + // Four labeled blocks, in the same order they are + // written: up/up, up/down Re, up/down Im, down/down. + std::vector uu(nm * nm); + std::vector re(nm * nm); + std::vector im(nm * nm); + std::vector dd(nm * nm); + + skip_compact_soc_label(ifdftu, word); + read_matrix_block(ifdftu, nm, uu.data()); + + skip_compact_soc_label(ifdftu, word); + read_matrix_block(ifdftu, nm, re.data()); + + // The imaginary block is needed by the PW (Pauli) + // path to reconstruct b2 = 2*Im(n_ud); the LCAO + // (real spin-basis) path discards it because its + // storage can not represent Im(n_ud). + skip_compact_soc_label(ifdftu, word); + read_matrix_block(ifdftu, nm, im.data()); + + skip_compact_soc_label(ifdftu, word); + read_matrix_block(ifdftu, nm, dd.data()); + + if (soc_layout == SOC_LAYOUT_PAULI) + { + // Reconstruct the 4 contiguous Pauli blocks + // [b0, b1, b2, b3] in the 2m x 2m flat buffer + // from the file's (n_uu, Re(n_ud), Im(n_ud), + // n_dd) representation. The writer prints + // n_uu = (b0 + b3)/2 + // Re(n_ud) = b1/2 + // Im(n_ud) = b2/2 + // n_dd = (b0 - b3)/2 + // so the inverse is + // b0 = n_uu + n_dd + // b1 = 2 * Re(n_ud) + // b2 = 2 * Im(n_ud) + // b3 = n_uu - n_dd. + ModuleBase::matrix& occ0 = occ.mat(iat, L, 0); + const int m2 = nm * nm; + for (int k = 0; k < m2; ++k) + { + occ0.c[0 * m2 + k] = uu[k] + dd[k]; + occ0.c[1 * m2 + k] = 2.0 * re[k]; + occ0.c[2 * m2 + k] = 2.0 * im[k]; + occ0.c[3 * m2 + k] = uu[k] - dd[k]; + } + } + else // SOC_LAYOUT_SPIN_BASIS_REAL + { + // Reconstruct the real 2m x 2m spin-basis matrix. + // Mirror the legacy writer: the Re block fills both + // (up, down) and (down, up) at the same (m0, m1). + // Im(n_ud) is dropped because the storage can not + // represent it. + for (int m0 = 0; m0 < nm; m0++) + { + for (int m1 = 0; m1 < nm; m1++) + { + const int k = m0 * nm + m1; + occ.set(iat, L, 0, m0, m1, uu[k]); + occ.set(iat, L, 0, m0, nm + m1, re[k]); + occ.set(iat, L, 0, nm + m0, m1, re[k]); + occ.set(iat, L, 0, nm + m0, nm + m1, dd[k]); + } + } + } + } + } + } + } + + return; +} + +} // namespace DFTU_BASE diff --git a/source/source_pw/module_pwdft/test/CMakeLists.txt b/source/source_pw/module_pwdft/test/CMakeLists.txt index ad7674a93e..a6cad027a3 100644 --- a/source/source_pw/module_pwdft/test/CMakeLists.txt +++ b/source/source_pw/module_pwdft/test/CMakeLists.txt @@ -28,6 +28,7 @@ AddTest( SOURCES dftu_base_test.cpp ../dftu_base.cpp ../dftu_base_io.cpp + ../dftu_base_io_read.cpp ../dftu_pw_tools.cpp ../yukawa_screening.cpp ../../../source_estate/occ_matrix.cpp From 052ee2bd42954657c0c2f8d1b07f85a859d0f54b Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Sat, 17 Oct 2026 20:35:22 +0800 Subject: [PATCH 02/13] Refactor DFT+U reader for strict compact input --- .../module_parameter/read_inp_exx_dftu.cpp | 2 +- source/source_pw/module_pwdft/dftu_base.cpp | 17 +- .../source_pw/module_pwdft/dftu_base_io.cpp | 1 - source/source_pw/module_pwdft/dftu_base_io.h | 6 +- .../module_pwdft/dftu_base_io_read.cpp | 620 +++++++++--------- source/source_pw/module_pwdft/dftu_types.h | 2 +- .../module_pwdft/test/dftu_base_test.cpp | 223 +++++++ 7 files changed, 560 insertions(+), 311 deletions(-) diff --git a/source/source_io/module_parameter/read_inp_exx_dftu.cpp b/source/source_io/module_parameter/read_inp_exx_dftu.cpp index bbecf5cace..5e4c416b92 100644 --- a/source/source_io/module_parameter/read_inp_exx_dftu.cpp +++ b/source/source_io/module_parameter/read_inp_exx_dftu.cpp @@ -816,7 +816,7 @@ void ReadInput::item_dftu() * 1: The first SCF step will use an initial density matrix read from a file, but for later steps, the onsite density matrix will be updated. * 2: The same onsite density matrix read from the file will be used throughout the entire calculation. -[NOTE] The occupation-matrix file is read from the directory specified by read_file_dir. The search order is dm_onsite_ini.txt -> occ_mat.txt -> dm_onsite.txt, so a previously output occ_mat.txt can be reused directly without renaming. The easiest way to obtain such a file is to run a DFT+U calculation with out_occ_mat=1 (the default); the resulting occ_mat.txt in the OUT.prefix directory (which always holds the occupation matrix of the latest electronic step) can be placed in read_file_dir as-is. When out_freq_ion is positive, separate occ_matg{geom_step}.txt files are also written, one per output ionic step; they may contain several electronic-step sections separated by comment headers. The format of the file is rather straight-forward.)"; +[NOTE] The occupation-matrix file is read from the directory specified by read_file_dir. The search order is dm_onsite_ini.txt -> occ_mat.txt. Input must use the current compact format written by ABACUS; legacy Atom=/L=/spin= files such as dm_onsite.txt are not supported. The easiest way to obtain a valid file is to run a DFT+U calculation with out_occ_mat=1 (the default); the resulting occ_mat.txt in the OUT.prefix directory (which always holds the occupation matrix of the latest electronic step) can be placed in read_file_dir as-is. When out_freq_ion is positive, separate occ_matg{geom_step}.txt files are also written, one per output ionic step; they may contain several electronic-step sections separated by comment headers.)"; item.default_value = "0"; item.unit = ""; read_sync_int(input.init_occ_mat); diff --git a/source/source_pw/module_pwdft/dftu_base.cpp b/source/source_pw/module_pwdft/dftu_base.cpp index 58caa93e78..f0dd2406e2 100644 --- a/source/source_pw/module_pwdft/dftu_base.cpp +++ b/source/source_pw/module_pwdft/dftu_base.cpp @@ -169,18 +169,17 @@ void Plus_U_Base::init_base(UnitCell& cell, else { // Try the user-prepared dm_onsite_ini.txt first, then fall back to - // the occupation-matrix snapshot written by a previous run - // (occ_mat.txt), and finally to the legacy dm_onsite.txt name. + // the compact occupation-matrix snapshot written by a previous + // run (occ_mat.txt). Legacy dm_onsite.txt input is not supported. const std::vector candidates = {"dm_onsite_ini.txt", - "occ_mat.txt", - "dm_onsite.txt"}; + "occ_mat.txt"}; const std::string readin_fn = DFTU_BASE::find_first_existing_file(global_readin_dir, candidates); if (readin_fn.empty()) { ModuleBase::WARNING_QUIT("Plus_U_Base::init_base", "init_occ_mat is set but no occupation-matrix file found in " + global_readin_dir - + ". Tried: dm_onsite_ini.txt, occ_mat.txt, dm_onsite.txt"); + + ". Tried: dm_onsite_ini.txt, occ_mat.txt"); } DFTU_BASE::read_occup_m(cell, this->occmat_, this->l_channel, this->init_occ_mat, readin_fn, init_chg, nspin, npol, soc_layout); @@ -202,10 +201,8 @@ void Plus_U_Base::init_base(UnitCell& cell, { if (init_chg == "file") { - // occ_mat.txt is the current output name; fall back to the - // legacy dm_onsite.txt so that output directories written by - // older versions can still be used for restarts. - const std::vector candidates = {"occ_mat.txt", "dm_onsite.txt"}; + // Only the current compact restart output is supported. + const std::vector candidates = {"occ_mat.txt"}; const std::string readin_fn = DFTU_BASE::find_first_existing_file(global_readin_dir, candidates); if (readin_fn.empty()) { @@ -213,7 +210,7 @@ void Plus_U_Base::init_base(UnitCell& cell, "init_chg is set to file but no occupation-matrix file " "found in " + global_readin_dir - + ". Tried: occ_mat.txt, dm_onsite.txt. Please do an scf " + + ". Tried: occ_mat.txt. Please do an scf " "calculation first."); } DFTU_BASE::read_occup_m(cell, this->occmat_, this->l_channel, this->init_occ_mat, diff --git a/source/source_pw/module_pwdft/dftu_base_io.cpp b/source/source_pw/module_pwdft/dftu_base_io.cpp index 05a13700ac..9449d8c5a9 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.cpp +++ b/source/source_pw/module_pwdft/dftu_base_io.cpp @@ -14,7 +14,6 @@ #include #include #include -#include #include #include #include diff --git a/source/source_pw/module_pwdft/dftu_base_io.h b/source/source_pw/module_pwdft/dftu_base_io.h index 8b24b39eed..179a283f93 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.h +++ b/source/source_pw/module_pwdft/dftu_base_io.h @@ -91,9 +91,9 @@ void append_ion_step_snapshot(const Plus_U_Base& dftu, /// Read the local occupation number matrix from file (rank 0 only). /// -/// The file format matches the output of write_occup_m(). When the file can -/// not be opened, the run quits with an error message that depends on -/// init_occ_mat and init_chg. +/// Only the current compact format written to occ_mat.txt is accepted. Legacy +/// Atom=/L=/spin= layouts and malformed or incomplete matrices terminate the +/// run with a diagnostic containing the filename and input line. /// /// @param soc_layout storage layout of the nspin == 4 occupation matrix; /// the PW path stores 4 contiguous Pauli blocks (b0, b1, b2, b3) and diff --git a/source/source_pw/module_pwdft/dftu_base_io_read.cpp b/source/source_pw/module_pwdft/dftu_base_io_read.cpp index a94d2ef2d8..5fc4ed47e5 100644 --- a/source/source_pw/module_pwdft/dftu_base_io_read.cpp +++ b/source/source_pw/module_pwdft/dftu_base_io_read.cpp @@ -4,47 +4,325 @@ #include "source_base/global_variable.h" #include "source_cell/unitcell.h" -#include +#include #include #include +#include +#include #include namespace { -// Consume the label of one SOC block in the compact layout, e.g. -// "spin 1 nelec 1.07922327" or "spin 12 re". The optional -// "nelec " token pair is skipped when present. -void skip_compact_soc_label(std::ifstream& ifs, char* word) +const char* const READ_CONTEXT = "DFTU_BASE::read_occup_m"; + +struct ReaderState +{ + ReaderState(std::ifstream& input, const std::string& filename) + : input(input), filename(filename), line_number(0) + { + } + + std::ifstream& input; + const std::string& filename; + int line_number; +}; + +[[noreturn]] void fail(const ReaderState& state, const std::string& message) +{ + std::ostringstream os; + os << "invalid compact occupation-matrix file " << state.filename; + if (state.line_number > 0) + { + os << " at line " << state.line_number; + } + os << ": " << message; + ModuleBase::WARNING_QUIT(READ_CONTEXT, os.str()); +} + +bool next_data_line(ReaderState& state, std::string& line) +{ + while (std::getline(state.input, line)) + { + ++state.line_number; + const std::string::size_type first = line.find_first_not_of(" \t\r"); + if (first == std::string::npos || line[first] == '#') + { + continue; + } + return true; + } + return false; +} + +std::string require_data_line(ReaderState& state, const std::string& expected) +{ + std::string line; + if (!next_data_line(state, line)) + { + fail(state, "unexpected end of file while reading " + expected); + } + return line; +} + +void require_line_end(ReaderState& state, + std::istringstream& line_stream, + const std::string& expected) +{ + std::string extra; + if (line_stream >> extra) + { + fail(state, "unexpected token '" + extra + "' after " + expected); + } +} + +double read_finite_double(ReaderState& state, + std::istringstream& line_stream, + const std::string& expected) +{ + double value = 0.0; + if (!(line_stream >> value) || !std::isfinite(value)) + { + fail(state, "expected a finite number for " + expected); + } + return value; +} + +void parse_atom_header(ReaderState& state, + const UnitCell& ucell, + const int expected_type, + const int expected_iat, + const int expected_l, + const int nspin) +{ + const std::string line = require_data_line(state, "an atom header"); + std::istringstream line_stream(line); + std::string element; + std::string atom_keyword; + std::string l_keyword; + std::string mag_keyword; + int atom_index = 0; + int l = -1; + if (!(line_stream >> element >> atom_keyword >> atom_index >> l_keyword >> l + >> mag_keyword)) + { + fail(state, "expected compact header ' Atom L mag ...'"); + } + if (atom_keyword != "Atom" || l_keyword != "L" || mag_keyword != "mag") + { + fail(state, "only the compact 'Atom L mag ...' header is supported"); + } + if (element != ucell.atoms[expected_type].label) + { + fail(state, "element label '" + element + "' does not match the unit cell"); + } + if (atom_index != expected_iat + 1) + { + fail(state, "atom index does not match the unit-cell order"); + } + if (l != expected_l) + { + fail(state, "L does not match the configured DFT+U channel"); + } + + const int magnetization_components = (nspin == 4) ? 3 : 1; + for (int component = 0; component < magnetization_components; ++component) + { + read_finite_double(state, line_stream, "magnetization"); + } + require_line_end(state, line_stream, "the compact atom header"); +} + +void parse_spin_nelec_label(ReaderState& state, const int expected_spin) +{ + const std::string line = require_data_line(state, "a spin block label"); + std::istringstream line_stream(line); + std::string spin_keyword; + std::string nelec_keyword; + int spin = 0; + if (!(line_stream >> spin_keyword >> spin >> nelec_keyword) + || spin_keyword != "spin" || spin != expected_spin + || nelec_keyword != "nelec") + { + fail(state, "expected compact label 'spin " + + std::to_string(expected_spin) + " nelec '"); + } + read_finite_double(state, line_stream, "nelec"); + require_line_end(state, line_stream, "the spin block label"); +} + +void parse_soc_offdiagonal_label(ReaderState& state, const std::string& component) +{ + const std::string line = require_data_line(state, "the spin 12 " + component + " label"); + std::istringstream line_stream(line); + std::string spin_keyword; + std::string block; + std::string actual_component; + if (!(line_stream >> spin_keyword >> block >> actual_component) + || spin_keyword != "spin" || block != "12" + || actual_component != component) + { + fail(state, "expected compact label 'spin 12 " + component + "'"); + } + require_line_end(state, line_stream, "the spin 12 " + component + " label"); +} + +std::vector parse_matrix(ReaderState& state, const int size) +{ + std::vector matrix(size * size, 0.0); + for (int row = 0; row < size; ++row) + { + const std::string line = require_data_line(state, "matrix row"); + std::istringstream line_stream(line); + for (int column = 0; column < size; ++column) + { + const int index = row * size + column; + matrix[index] = read_finite_double(state, line_stream, "matrix element"); + } + require_line_end(state, line_stream, "a matrix row"); + } + return matrix; +} + +void read_collinear_atom(ReaderState& state, + OccupationMatrix& occ, + const int iat, + const int l) +{ + const int size = 2 * l + 1; + std::vector> spin_matrices(2); + for (int spin = 0; spin < 2; ++spin) + { + parse_spin_nelec_label(state, spin + 1); + spin_matrices[spin] = parse_matrix(state, size); + } + + for (int spin = 0; spin < 2; ++spin) + { + for (int row = 0; row < size; ++row) + { + for (int column = 0; column < size; ++column) + { + const int index = row * size + column; + occ.set(iat, l, spin, row, column, spin_matrices[spin][index]); + } + } + } +} + +void store_soc_pauli(OccupationMatrix& occ, + const int iat, + const int l, + const std::vector>& blocks) { - ifs >> word; // "spin" - ifs >> word; // block index: "1", "2" or "12" - ifs >> word; // "nelec", "re" or "im" - if (strcmp(word, "nelec") == 0) + ModuleBase::matrix& matrix = occ.mat(iat, l, 0); + const int block_size = static_cast(blocks[0].size()); + for (int index = 0; index < block_size; ++index) { - double nelec = 0.0; - ifs >> nelec; + matrix.c[0 * block_size + index] = blocks[0][index] + blocks[3][index]; + matrix.c[1 * block_size + index] = 2.0 * blocks[1][index]; + matrix.c[2 * block_size + index] = 2.0 * blocks[2][index]; + matrix.c[3 * block_size + index] = blocks[0][index] - blocks[3][index]; } } -// Read n*n doubles of one compact matrix block row by row. When out is -// non-null the values are stored, otherwise they are read and discarded. -void read_matrix_block(std::ifstream& ifs, const int n, double* out) +void store_soc_spin_basis(OccupationMatrix& occ, + const int iat, + const int l, + const std::vector>& blocks) { - for (int i = 0; i < n * n; ++i) + const int size = 2 * l + 1; + for (int row = 0; row < size; ++row) { - double value = 0.0; - ifs >> value; - if (ifs.fail()) + for (int column = 0; column < size; ++column) { - ModuleBase::WARNING_QUIT("read_matrix_block", - "failed to read a matrix element from the occupation matrix file"); + const int index = row * size + column; + occ.set(iat, l, 0, row, column, blocks[0][index]); + occ.set(iat, l, 0, row, size + column, blocks[1][index]); + occ.set(iat, l, 0, size + row, column, blocks[1][index]); + occ.set(iat, l, 0, size + row, size + column, blocks[3][index]); } - if (out != nullptr) + } +} + +void read_soc_atom(ReaderState& state, + OccupationMatrix& occ, + const int iat, + const int l, + const DFTU_BASE::OccmatSocLayout soc_layout) +{ + const int size = 2 * l + 1; + std::vector> blocks(4); + parse_spin_nelec_label(state, 1); + blocks[0] = parse_matrix(state, size); + parse_soc_offdiagonal_label(state, "re"); + blocks[1] = parse_matrix(state, size); + parse_soc_offdiagonal_label(state, "im"); + blocks[2] = parse_matrix(state, size); + parse_spin_nelec_label(state, 2); + blocks[3] = parse_matrix(state, size); + + if (soc_layout == DFTU_BASE::SOC_LAYOUT_PAULI) + { + store_soc_pauli(occ, iat, l, blocks); + } + else + { + store_soc_spin_basis(occ, iat, l, blocks); + } +} + +void validate_reader_configuration(const ReaderState& state, + const UnitCell& ucell, + const std::vector& l_channel, + const int nspin, + const int npol, + const DFTU_BASE::OccmatSocLayout soc_layout) +{ + if (static_cast(l_channel.size()) != ucell.ntype) + { + fail(state, "l_channel size does not match the number of atom types"); + } + for (int type = 0; type < ucell.ntype; ++type) + { + if (l_channel[type] < -1 || l_channel[type] > ucell.atoms[type].nwl) { - out[i] = value; + fail(state, "a configured DFT+U channel is outside the unit-cell range"); } } + if (nspin != 1 && nspin != 2 && nspin != 4) + { + fail(state, "nspin must be 1, 2 or 4"); + } + const int expected_npol = (nspin == 4) ? 2 : 1; + if (npol != expected_npol) + { + fail(state, "npol is inconsistent with nspin"); + } + if (soc_layout != DFTU_BASE::SOC_LAYOUT_PAULI + && soc_layout != DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL) + { + fail(state, "unknown SOC occupation-matrix storage layout"); + } +} + +void report_open_error(const int init_occ_mat, + const std::string& filename, + const std::string& init_chg) +{ + if (init_occ_mat > 0) + { + ModuleBase::WARNING_QUIT(READ_CONTEXT, + "Can not open compact occupation-matrix file " + filename); + } + if (init_chg == "file") + { + ModuleBase::WARNING_QUIT(READ_CONTEXT, + "Can not find the compact occupation-matrix file " + + filename + ". Please do scf calculation first"); + } + ModuleBase::WARNING_QUIT(READ_CONTEXT, + "Can not open compact occupation-matrix file " + filename); } } // namespace @@ -63,305 +341,57 @@ void read_occup_m(const UnitCell& ucell, OccmatSocLayout soc_layout) { ModuleBase::TITLE("DFTU_BASE", "read_occup_m"); - if (GlobalV::MY_RANK != 0) { return; } - std::ifstream ifdftu(fn.c_str(), std::ios::in); - - if (!ifdftu) + std::ifstream input(fn.c_str()); + if (!input) { - if (init_occ_mat > 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "Can not find the file dm_onsite_ini.txt. Please check your dm_onsite_ini.txt"); - } - else - { - if (init_chg == "file") - { - const std::string not_found_msg = "Can not find the file " + fn - + ". Please do scf calculation first"; - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", not_found_msg); - } - } - const std::string open_failed_msg = "Can not open " + fn + " file"; - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", open_failed_msg); + report_open_error(init_occ_mat, fn, init_chg); } GlobalV::ofs_running << " DFT+U: read occupation matrix from " << fn << std::endl; std::cout << " DFT+U: read occupation matrix from " << fn << std::endl; - ifdftu.clear(); - ifdftu.seekg(0); - - char word[20]; + ReaderState state(input, fn); + validate_reader_configuration(state, ucell, l_channel, nspin, npol, soc_layout); - int T = 0; - int iat = 0; - int spin = 0; - int L = 0; - int zeta = 0; - - while (true) + int atom_count = 0; + for (int type = 0; type < ucell.ntype; ++type) { - if (!(ifdftu >> word)) - { - break; - } - - // Comment lines ("# ...") carry the provenance header and the - // electronic-step marker in the compact layout; skip the line. - if (word[0] == '#') - { - ifdftu.ignore(150, '\n'); - continue; - } - - // The legacy "Atom=" token and the compact "Atom" label both - // start one atom block; every other token (element prefix, log - // sentences) is ignored. - const bool legacy_atom = (strcmp(word, "Atom=") == 0); - const bool compact_atom = (strcmp(word, "Atom") == 0); - if (!legacy_atom && !compact_atom) + const int l = l_channel[type]; + if (l < 0) { continue; } - const bool legacy_layout = legacy_atom; - - ifdftu >> iat; - iat -= 1; - - if (legacy_layout) + for (int atom = 0; atom < ucell.atoms[type].na; ++atom) { - ifdftu >> word; - - // Accept both the compact "L=2" and the legacy split "L= 2". - if (strncmp("L=", word, 2) != 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); - } - if (strlen(word) > 2) - { - L = atoi(word + 2); - } - else - { - ifdftu >> L; - } - ifdftu >> word; - - // The ORBITAL= token is optional in newly written files; when - // absent the radial channel defaults to 0. - if (strncmp("ORBITAL=", word, 8) == 0) + const int iat = ucell.itia2iat(type, atom); + parse_atom_header(state, ucell, type, iat, l, nspin); + if (nspin == 4) { - if (strlen(word) > 8) - { - zeta = atoi(word + 8); - } - else - { - ifdftu >> zeta; - } + read_soc_atom(state, occ, iat, l, soc_layout); } else { - zeta = 0; - } - ifdftu.ignore(150, '\n'); - - if (zeta != 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", - "only the first radial channel (ORBITAL=0) is supported"); + read_collinear_atom(state, occ, iat, l); } + ++atom_count; } - else - { - // Compact header, e.g. " Fe Atom 1 L 2 mag 0.03..." - ifdftu >> word; - if (strcmp(word, "L") != 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); - } - ifdftu >> L; - // The optional "mag ..." tail is redundant for a restart. - ifdftu.ignore(150, '\n'); - zeta = 0; - } - - T = ucell.iat2it[iat]; - const int NL = ucell.atoms[T].nwl + 1; - - for (int l = 0; l < NL; l++) - { - if (l != l_channel[T]) - { - continue; - } - - const int nm = 2 * L + 1; + } - if (nspin == 1 || nspin == 2) - { - for (int is = 0; is < 2; is++) - { - if (legacy_layout) - { - ifdftu >> word; - if (strcmp("spin=", word) != 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); - } - ifdftu >> spin; - spin -= 1; - ifdftu.ignore(150, '\n'); - - double value = 0.0; - for (int m0 = 0; m0 < nm; m0++) - { - for (int m1 = 0; m1 < nm; m1++) - { - ifdftu >> value; - occ.set(iat, L, spin, m0, m1, value); - } - ifdftu.ignore(150, '\n'); - } - } - else - { - // Compact label: "spin 1 nelec 0.468..." - ifdftu >> word; - if (strcmp("spin", word) != 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); - } - int spin_idx = 0; - ifdftu >> spin_idx; - if (spin_idx != is + 1) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG SPIN INDEX IN Plus_U FILE"); - } - ifdftu >> word; - if (strcmp("nelec", word) != 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); - } - double nelec = 0.0; - ifdftu >> nelec; - - double value = 0.0; - for (int m0 = 0; m0 < nm; m0++) - { - for (int m1 = 0; m1 < nm; m1++) - { - ifdftu >> value; - occ.set(iat, L, is, m0, m1, value); - } - } - } - } - } - else if (nspin == 4) // SOC - { - if (legacy_layout) - { - double value = 0.0; - for (int m0 = 0; m0 < nm; m0++) - { - for (int ipol0 = 0; ipol0 < npol; ipol0++) - { - const int m0_all = m0 + nm * ipol0; - - for (int m1 = 0; m1 < nm; m1++) - { - for (int ipol1 = 0; ipol1 < npol; ipol1++) - { - int m1_all = m1 + nm * ipol1; - ifdftu >> value; - occ.set(iat, L, 0, m0_all, m1_all, value); - } - } - ifdftu.ignore(150, '\n'); - } - } - } - else - { - // Four labeled blocks, in the same order they are - // written: up/up, up/down Re, up/down Im, down/down. - std::vector uu(nm * nm); - std::vector re(nm * nm); - std::vector im(nm * nm); - std::vector dd(nm * nm); - - skip_compact_soc_label(ifdftu, word); - read_matrix_block(ifdftu, nm, uu.data()); - - skip_compact_soc_label(ifdftu, word); - read_matrix_block(ifdftu, nm, re.data()); - - // The imaginary block is needed by the PW (Pauli) - // path to reconstruct b2 = 2*Im(n_ud); the LCAO - // (real spin-basis) path discards it because its - // storage can not represent Im(n_ud). - skip_compact_soc_label(ifdftu, word); - read_matrix_block(ifdftu, nm, im.data()); - - skip_compact_soc_label(ifdftu, word); - read_matrix_block(ifdftu, nm, dd.data()); - - if (soc_layout == SOC_LAYOUT_PAULI) - { - // Reconstruct the 4 contiguous Pauli blocks - // [b0, b1, b2, b3] in the 2m x 2m flat buffer - // from the file's (n_uu, Re(n_ud), Im(n_ud), - // n_dd) representation. The writer prints - // n_uu = (b0 + b3)/2 - // Re(n_ud) = b1/2 - // Im(n_ud) = b2/2 - // n_dd = (b0 - b3)/2 - // so the inverse is - // b0 = n_uu + n_dd - // b1 = 2 * Re(n_ud) - // b2 = 2 * Im(n_ud) - // b3 = n_uu - n_dd. - ModuleBase::matrix& occ0 = occ.mat(iat, L, 0); - const int m2 = nm * nm; - for (int k = 0; k < m2; ++k) - { - occ0.c[0 * m2 + k] = uu[k] + dd[k]; - occ0.c[1 * m2 + k] = 2.0 * re[k]; - occ0.c[2 * m2 + k] = 2.0 * im[k]; - occ0.c[3 * m2 + k] = uu[k] - dd[k]; - } - } - else // SOC_LAYOUT_SPIN_BASIS_REAL - { - // Reconstruct the real 2m x 2m spin-basis matrix. - // Mirror the legacy writer: the Re block fills both - // (up, down) and (down, up) at the same (m0, m1). - // Im(n_ud) is dropped because the storage can not - // represent it. - for (int m0 = 0; m0 < nm; m0++) - { - for (int m1 = 0; m1 < nm; m1++) - { - const int k = m0 * nm + m1; - occ.set(iat, L, 0, m0, m1, uu[k]); - occ.set(iat, L, 0, m0, nm + m1, re[k]); - occ.set(iat, L, 0, nm + m0, m1, re[k]); - occ.set(iat, L, 0, nm + m0, nm + m1, dd[k]); - } - } - } - } - } - } + if (atom_count == 0) + { + fail(state, "the unit cell has no configured DFT+U atom blocks"); } - return; + std::string extra_line; + if (next_data_line(state, extra_line)) + { + fail(state, "unexpected content after the final atom block"); + } } } // namespace DFTU_BASE diff --git a/source/source_pw/module_pwdft/dftu_types.h b/source/source_pw/module_pwdft/dftu_types.h index f409e5b884..e1f1bee13d 100644 --- a/source/source_pw/module_pwdft/dftu_types.h +++ b/source/source_pw/module_pwdft/dftu_types.h @@ -26,7 +26,7 @@ struct OccmatOutputCfg /// Text format used by write_occup_m(). enum OccmatTextFormat { - OCMAT_FMT_LEGACY, ///< whitespace-token layout parsed by read_occup_m() + OCMAT_FMT_LEGACY, ///< legacy diagnostic layout; not accepted on input OCMAT_FMT_READABLE ///< compact human-readable layout for snapshot files }; 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 aca4b5f0e4..66936fabde 100644 --- a/source/source_pw/module_pwdft/test/dftu_base_test.cpp +++ b/source/source_pw/module_pwdft/test/dftu_base_test.cpp @@ -107,6 +107,38 @@ class DFTUBaseTest : public testing::Test } }; +namespace +{ + +void write_constant_matrix(std::ostream& os, + const int size, + const double value, + const int first_row_columns) +{ + for (int row = 0; row < size; ++row) + { + const int columns = (row == 0 && first_row_columns >= 0) + ? first_row_columns + : size; + for (int column = 0; column < columns; ++column) + { + os << " " << value; + } + os << "\n"; + } +} + +void write_compact_collinear_body(std::ostream& os) +{ + os << "Fe Atom 1 L 2 mag 0.0\n"; + os << "spin 1 nelec 0.5\n"; + write_constant_matrix(os, 5, 0.1, 5); + os << "spin 2 nelec 0.5\n"; + write_constant_matrix(os, 5, 0.2, 5); +} + +} // namespace + /// After a true -> false re-initialization the Yukawa object must be /// released so that use_yukawa() reflects the latest argument. TEST_F(DFTUBaseTest, InitBaseYukawaTrueThenFalseClearsState) @@ -780,6 +812,197 @@ TEST_F(DFTUBaseTest, InitBaseReadsOccMatSocPauliRoundtrip) } } +TEST_F(DFTUBaseTest, ReadOccMatCompactHandlesLongComment) +{ + const std::string fn = testing::TempDir() + "compact_long_comment.txt"; + { + std::ofstream ofs(fn); + ASSERT_TRUE(ofs.is_open()); + ofs << "# " << std::string(300, 'x') << "\n"; + write_compact_collinear_body(ofs); + } + + const std::vector l_channel = {2}; + OccupationMatrix occ; + occ.init(ucell, l_channel, 2, 1); + DFTU_BASE::read_occup_m(ucell, + occ, + l_channel, + 2, + fn, + "none", + 2, + 1, + DFTU_BASE::SOC_LAYOUT_PAULI); + EXPECT_DOUBLE_EQ(occ.get(0, 2, 0, 0, 0), 0.1); + EXPECT_DOUBLE_EQ(occ.get(0, 2, 1, 4, 4), 0.2); +} + +TEST_F(DFTUBaseTest, ReadOccMatCompactSocRealSpinBasis) +{ + const std::string fn = testing::TempDir() + "compact_soc_real.txt"; + { + std::ofstream ofs(fn); + ASSERT_TRUE(ofs.is_open()); + ofs << "Fe Atom 1 L 2 mag 0.0 0.0 0.0\n"; + ofs << "spin 1 nelec 0.5\n"; + write_constant_matrix(ofs, 5, 1.0, 5); + ofs << "spin 12 re\n"; + write_constant_matrix(ofs, 5, 2.0, 5); + ofs << "spin 12 im\n"; + write_constant_matrix(ofs, 5, 3.0, 5); + ofs << "spin 2 nelec 0.5\n"; + write_constant_matrix(ofs, 5, 4.0, 5); + } + + const std::vector l_channel = {2}; + OccupationMatrix occ; + occ.init(ucell, l_channel, 4, 2); + DFTU_BASE::read_occup_m(ucell, + occ, + l_channel, + 2, + fn, + "none", + 4, + 2, + DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); + EXPECT_DOUBLE_EQ(occ.get(0, 2, 0, 0, 0), 1.0); + EXPECT_DOUBLE_EQ(occ.get(0, 2, 0, 0, 5), 2.0); + EXPECT_DOUBLE_EQ(occ.get(0, 2, 0, 5, 0), 2.0); + EXPECT_DOUBLE_EQ(occ.get(0, 2, 0, 5, 5), 4.0); +} + +TEST_F(DFTUBaseTest, ReadOccMatRejectsLegacyHeader) +{ + const std::string fn = testing::TempDir() + "legacy_header.txt"; + { + std::ofstream ofs(fn); + ASSERT_TRUE(ofs.is_open()); + ofs << "Fe 1 Atom=1 L=2 0 0 0\n"; + } + const std::vector l_channel = {2}; + OccupationMatrix occ; + occ.init(ucell, l_channel, 2, 1); + EXPECT_EXIT(DFTU_BASE::read_occup_m(ucell, + occ, + l_channel, + 2, + fn, + "none", + 2, + 1, + DFTU_BASE::SOC_LAYOUT_PAULI), + testing::ExitedWithCode(1), + ""); +} + +TEST_F(DFTUBaseTest, ReadOccMatRejectsShortMatrixRow) +{ + const std::string fn = testing::TempDir() + "short_matrix_row.txt"; + { + std::ofstream ofs(fn); + ASSERT_TRUE(ofs.is_open()); + ofs << "Fe Atom 1 L 2 mag 0.0\n"; + ofs << "spin 1 nelec 0.5\n"; + write_constant_matrix(ofs, 5, 0.1, 4); + } + const std::vector l_channel = {2}; + OccupationMatrix occ; + occ.init(ucell, l_channel, 2, 1); + EXPECT_EXIT(DFTU_BASE::read_occup_m(ucell, + occ, + l_channel, + 2, + fn, + "none", + 2, + 1, + DFTU_BASE::SOC_LAYOUT_PAULI), + testing::ExitedWithCode(1), + ""); +} + +TEST_F(DFTUBaseTest, ReadOccMatRejectsLongMatrixRow) +{ + const std::string fn = testing::TempDir() + "long_matrix_row.txt"; + { + std::ofstream ofs(fn); + ASSERT_TRUE(ofs.is_open()); + ofs << "Fe Atom 1 L 2 mag 0.0\n"; + ofs << "spin 1 nelec 0.5\n"; + write_constant_matrix(ofs, 5, 0.1, 6); + } + const std::vector l_channel = {2}; + OccupationMatrix occ; + occ.init(ucell, l_channel, 2, 1); + EXPECT_EXIT(DFTU_BASE::read_occup_m(ucell, + occ, + l_channel, + 2, + fn, + "none", + 2, + 1, + DFTU_BASE::SOC_LAYOUT_PAULI), + testing::ExitedWithCode(1), + ""); +} + +TEST_F(DFTUBaseTest, ReadOccMatRejectsMissingSocBlock) +{ + const std::string fn = testing::TempDir() + "missing_soc_block.txt"; + { + std::ofstream ofs(fn); + ASSERT_TRUE(ofs.is_open()); + ofs << "Fe Atom 1 L 2 mag 0.0 0.0 0.0\n"; + ofs << "spin 1 nelec 0.5\n"; + write_constant_matrix(ofs, 5, 1.0, 5); + ofs << "spin 12 re\n"; + write_constant_matrix(ofs, 5, 2.0, 5); + ofs << "spin 2 nelec 0.5\n"; + } + const std::vector l_channel = {2}; + OccupationMatrix occ; + occ.init(ucell, l_channel, 4, 2); + EXPECT_EXIT(DFTU_BASE::read_occup_m(ucell, + occ, + l_channel, + 2, + fn, + "none", + 4, + 2, + DFTU_BASE::SOC_LAYOUT_PAULI), + testing::ExitedWithCode(1), + ""); +} + +TEST_F(DFTUBaseTest, ReadOccMatRejectsExtraContent) +{ + const std::string fn = testing::TempDir() + "extra_content.txt"; + { + std::ofstream ofs(fn); + ASSERT_TRUE(ofs.is_open()); + write_compact_collinear_body(ofs); + ofs << "unexpected_trailing_token_that_is_longer_than_twenty_characters\n"; + } + const std::vector l_channel = {2}; + OccupationMatrix occ; + occ.init(ucell, l_channel, 2, 1); + EXPECT_EXIT(DFTU_BASE::read_occup_m(ucell, + occ, + l_channel, + 2, + fn, + "none", + 2, + 1, + DFTU_BASE::SOC_LAYOUT_PAULI), + testing::ExitedWithCode(1), + ""); +} + /// For init_occ_mat=1 the file-loaded occupation matrix must seed the /// Hamiltonian at the first SCF step. Before the fix the whole /// iter_init_dftu_pw was skipped at the first step, so pot_onsite and From 04634f54e63a7d8d4bf5cb34287239b2a9866f3e Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Sat, 17 Oct 2026 20:41:09 +0800 Subject: [PATCH 03/13] Rename and simplify DFT+U compact reader --- .../module_dftu/unittests/CMakeLists.txt | 4 +-- source/source_pw/module_pwdft/CMakeLists.txt | 2 +- source/source_pw/module_pwdft/dftu_base.cpp | 8 +++--- .../source_pw/module_pwdft/dftu_base_io.cpp | 13 +++++----- source/source_pw/module_pwdft/dftu_base_io.h | 2 -- ...tu_base_io_read.cpp => dftu_base_read.cpp} | 25 ++++--------------- .../module_pwdft/test/CMakeLists.txt | 2 +- .../module_pwdft/test/dftu_base_test.cpp | 14 ----------- 8 files changed, 19 insertions(+), 51 deletions(-) rename source/source_pw/module_pwdft/{dftu_base_io_read.cpp => dftu_base_read.cpp} (92%) diff --git a/source/source_lcao/module_dftu/unittests/CMakeLists.txt b/source/source_lcao/module_dftu/unittests/CMakeLists.txt index 6a50a09ffe..56d98da814 100644 --- a/source/source_lcao/module_dftu/unittests/CMakeLists.txt +++ b/source/source_lcao/module_dftu/unittests/CMakeLists.txt @@ -91,7 +91,7 @@ AddTest( ../../../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_read.cpp ../../../source_pw/module_pwdft/yukawa_screening.cpp ../../../source_estate/occ_matrix.cpp ../../../source_estate/occ_mixer.cpp @@ -137,7 +137,7 @@ AddTest( ../../../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_read.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 404fb05b4f..0e07989054 100644 --- a/source/source_pw/module_pwdft/CMakeLists.txt +++ b/source/source_pw/module_pwdft/CMakeLists.txt @@ -31,7 +31,7 @@ endif() list(APPEND objects dftu_base.cpp dftu_base_io.cpp - dftu_base_io_read.cpp + dftu_base_read.cpp dftu_pw.cpp dftu_pw_tools.cpp yukawa_screening.cpp diff --git a/source/source_pw/module_pwdft/dftu_base.cpp b/source/source_pw/module_pwdft/dftu_base.cpp index f0dd2406e2..0bb70b999a 100644 --- a/source/source_pw/module_pwdft/dftu_base.cpp +++ b/source/source_pw/module_pwdft/dftu_base.cpp @@ -181,8 +181,8 @@ void Plus_U_Base::init_base(UnitCell& cell, + global_readin_dir + ". Tried: dm_onsite_ini.txt, occ_mat.txt"); } - DFTU_BASE::read_occup_m(cell, this->occmat_, this->l_channel, this->init_occ_mat, - readin_fn, init_chg, nspin, npol, soc_layout); + 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); #endif @@ -213,8 +213,8 @@ void Plus_U_Base::init_base(UnitCell& cell, + ". Tried: occ_mat.txt. Please do an scf " "calculation first."); } - DFTU_BASE::read_occup_m(cell, this->occmat_, this->l_channel, this->init_occ_mat, - readin_fn, init_chg, nspin, npol, soc_layout); + 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); #endif diff --git a/source/source_pw/module_pwdft/dftu_base_io.cpp b/source/source_pw/module_pwdft/dftu_base_io.cpp index 9449d8c5a9..a8f680b33e 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.cpp +++ b/source/source_pw/module_pwdft/dftu_base_io.cpp @@ -425,8 +425,8 @@ void write_snapshot_section(std::ostream& os, // 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 - // zero matrix. read_occup_m() skips the token like any other non-atom - // word, so the file stays parseable on restart. + // zero matrix. Such an incomplete snapshot is diagnostic output and is + // intentionally rejected by the strict restart reader. if (!occmat_ready) { os << "\n N/A\n"; @@ -637,9 +637,8 @@ void write_occup_m(const Plus_U_Base& dftu, } else { - // Legacy token layout parsed by read_occup_m(); the - // " " prefix is skipped when scanning - // for "Atom=". + // 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; @@ -951,8 +950,8 @@ void write_occup_m(const Plus_U_Base& dftu, } else { - // Real 2m x 2m matrix in spin basis; this is the - // layout parsed by read_occup_m() on restart. + // 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) diff --git a/source/source_pw/module_pwdft/dftu_base_io.h b/source/source_pw/module_pwdft/dftu_base_io.h index 179a283f93..babcc76908 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.h +++ b/source/source_pw/module_pwdft/dftu_base_io.h @@ -103,9 +103,7 @@ void append_ion_step_snapshot(const Plus_U_Base& dftu, void read_occup_m(const UnitCell& ucell, OccupationMatrix& occ, const std::vector& l_channel, - const int init_occ_mat, const std::string& fn, - const std::string& init_chg, int nspin, int npol, OccmatSocLayout soc_layout); diff --git a/source/source_pw/module_pwdft/dftu_base_io_read.cpp b/source/source_pw/module_pwdft/dftu_base_read.cpp similarity index 92% rename from source/source_pw/module_pwdft/dftu_base_io_read.cpp rename to source/source_pw/module_pwdft/dftu_base_read.cpp index 5fc4ed47e5..c887b260be 100644 --- a/source/source_pw/module_pwdft/dftu_base_io_read.cpp +++ b/source/source_pw/module_pwdft/dftu_base_read.cpp @@ -14,7 +14,7 @@ namespace { -const char* const READ_CONTEXT = "DFTU_BASE::read_occup_m"; +const char* const read_context = "DFTU_BASE::read_occup_m"; struct ReaderState { @@ -37,7 +37,7 @@ struct ReaderState os << " at line " << state.line_number; } os << ": " << message; - ModuleBase::WARNING_QUIT(READ_CONTEXT, os.str()); + ModuleBase::WARNING_QUIT(read_context, os.str()); } bool next_data_line(ReaderState& state, std::string& line) @@ -306,22 +306,9 @@ void validate_reader_configuration(const ReaderState& state, } } -void report_open_error(const int init_occ_mat, - const std::string& filename, - const std::string& init_chg) +void report_open_error(const std::string& filename) { - if (init_occ_mat > 0) - { - ModuleBase::WARNING_QUIT(READ_CONTEXT, - "Can not open compact occupation-matrix file " + filename); - } - if (init_chg == "file") - { - ModuleBase::WARNING_QUIT(READ_CONTEXT, - "Can not find the compact occupation-matrix file " - + filename + ". Please do scf calculation first"); - } - ModuleBase::WARNING_QUIT(READ_CONTEXT, + ModuleBase::WARNING_QUIT(read_context, "Can not open compact occupation-matrix file " + filename); } @@ -333,9 +320,7 @@ namespace DFTU_BASE void read_occup_m(const UnitCell& ucell, OccupationMatrix& occ, const std::vector& l_channel, - const int init_occ_mat, const std::string& fn, - const std::string& init_chg, int nspin, int npol, OccmatSocLayout soc_layout) @@ -349,7 +334,7 @@ void read_occup_m(const UnitCell& ucell, std::ifstream input(fn.c_str()); if (!input) { - report_open_error(init_occ_mat, fn, init_chg); + report_open_error(fn); } GlobalV::ofs_running << " DFT+U: read occupation matrix from " << fn << std::endl; diff --git a/source/source_pw/module_pwdft/test/CMakeLists.txt b/source/source_pw/module_pwdft/test/CMakeLists.txt index a6cad027a3..51ebf66857 100644 --- a/source/source_pw/module_pwdft/test/CMakeLists.txt +++ b/source/source_pw/module_pwdft/test/CMakeLists.txt @@ -28,7 +28,7 @@ AddTest( SOURCES dftu_base_test.cpp ../dftu_base.cpp ../dftu_base_io.cpp - ../dftu_base_io_read.cpp + ../dftu_base_read.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 66936fabde..3d96a1137d 100644 --- a/source/source_pw/module_pwdft/test/dftu_base_test.cpp +++ b/source/source_pw/module_pwdft/test/dftu_base_test.cpp @@ -828,9 +828,7 @@ TEST_F(DFTUBaseTest, ReadOccMatCompactHandlesLongComment) DFTU_BASE::read_occup_m(ucell, occ, l_channel, - 2, fn, - "none", 2, 1, DFTU_BASE::SOC_LAYOUT_PAULI); @@ -861,9 +859,7 @@ TEST_F(DFTUBaseTest, ReadOccMatCompactSocRealSpinBasis) DFTU_BASE::read_occup_m(ucell, occ, l_channel, - 2, fn, - "none", 4, 2, DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); @@ -887,9 +883,7 @@ TEST_F(DFTUBaseTest, ReadOccMatRejectsLegacyHeader) EXPECT_EXIT(DFTU_BASE::read_occup_m(ucell, occ, l_channel, - 2, fn, - "none", 2, 1, DFTU_BASE::SOC_LAYOUT_PAULI), @@ -913,9 +907,7 @@ TEST_F(DFTUBaseTest, ReadOccMatRejectsShortMatrixRow) EXPECT_EXIT(DFTU_BASE::read_occup_m(ucell, occ, l_channel, - 2, fn, - "none", 2, 1, DFTU_BASE::SOC_LAYOUT_PAULI), @@ -939,9 +931,7 @@ TEST_F(DFTUBaseTest, ReadOccMatRejectsLongMatrixRow) EXPECT_EXIT(DFTU_BASE::read_occup_m(ucell, occ, l_channel, - 2, fn, - "none", 2, 1, DFTU_BASE::SOC_LAYOUT_PAULI), @@ -968,9 +958,7 @@ TEST_F(DFTUBaseTest, ReadOccMatRejectsMissingSocBlock) EXPECT_EXIT(DFTU_BASE::read_occup_m(ucell, occ, l_channel, - 2, fn, - "none", 4, 2, DFTU_BASE::SOC_LAYOUT_PAULI), @@ -993,9 +981,7 @@ TEST_F(DFTUBaseTest, ReadOccMatRejectsExtraContent) EXPECT_EXIT(DFTU_BASE::read_occup_m(ucell, occ, l_channel, - 2, fn, - "none", 2, 1, DFTU_BASE::SOC_LAYOUT_PAULI), From 686c197b27fe95106fc284ff8535631446a0d859 Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Sat, 17 Oct 2026 21:21:52 +0800 Subject: [PATCH 04/13] Fix legacy Makefile linkage for DFT+U reader --- source/Makefile.Objects | 1 + 1 file changed, 1 insertion(+) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 4956896a05..6083ff86fa 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -951,6 +951,7 @@ OBJS_SRCPW=h_ewald_pw.o\ update_cell_pw.o\ dftu_base.o\ dftu_base_io.o\ + dftu_base_read.o\ dftu_pw.o\ dftu_pw_tools.o\ yukawa_screening.o\ From 9b45b0e1d4893695313590fe49b454f03c0dd41b Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Sat, 17 Oct 2026 21:22:44 +0800 Subject: [PATCH 05/13] Scope DFT+U occupation reader helpers explicitly --- .../source_pw/module_pwdft/dftu_base_read.cpp | 216 +++++++++--------- 1 file changed, 110 insertions(+), 106 deletions(-) diff --git a/source/source_pw/module_pwdft/dftu_base_read.cpp b/source/source_pw/module_pwdft/dftu_base_read.cpp index c887b260be..1874f6ebe9 100644 --- a/source/source_pw/module_pwdft/dftu_base_read.cpp +++ b/source/source_pw/module_pwdft/dftu_base_read.cpp @@ -11,14 +11,14 @@ #include #include -namespace +namespace DFTU_BASE { -const char* const read_context = "DFTU_BASE::read_occup_m"; +const char* const read_occup_context = "DFTU_BASE::read_occup_m"; -struct ReaderState +struct OccMatReadState { - ReaderState(std::ifstream& input, const std::string& filename) + OccMatReadState(std::ifstream& input, const std::string& filename) : input(input), filename(filename), line_number(0) { } @@ -28,7 +28,8 @@ struct ReaderState int line_number; }; -[[noreturn]] void fail(const ReaderState& state, const std::string& message) +[[noreturn]] static void read_occup_fail(const OccMatReadState& state, + const std::string& message) { std::ostringstream os; os << "invalid compact occupation-matrix file " << state.filename; @@ -37,10 +38,10 @@ struct ReaderState os << " at line " << state.line_number; } os << ": " << message; - ModuleBase::WARNING_QUIT(read_context, os.str()); + ModuleBase::WARNING_QUIT(read_occup_context, os.str()); } -bool next_data_line(ReaderState& state, std::string& line) +static bool read_occup_next_data_line(OccMatReadState& state, std::string& line) { while (std::getline(state.input, line)) { @@ -55,47 +56,48 @@ bool next_data_line(ReaderState& state, std::string& line) return false; } -std::string require_data_line(ReaderState& state, const std::string& expected) +static std::string read_occup_require_data_line(OccMatReadState& state, + const std::string& expected) { std::string line; - if (!next_data_line(state, line)) + if (!read_occup_next_data_line(state, line)) { - fail(state, "unexpected end of file while reading " + expected); + read_occup_fail(state, "unexpected end of file while reading " + expected); } return line; } -void require_line_end(ReaderState& state, - std::istringstream& line_stream, - const std::string& expected) +static void read_occup_require_line_end(OccMatReadState& state, + std::istringstream& line_stream, + const std::string& expected) { std::string extra; if (line_stream >> extra) { - fail(state, "unexpected token '" + extra + "' after " + expected); + read_occup_fail(state, "unexpected token '" + extra + "' after " + expected); } } -double read_finite_double(ReaderState& state, - std::istringstream& line_stream, - const std::string& expected) +static double read_occup_finite_double(OccMatReadState& state, + std::istringstream& line_stream, + const std::string& expected) { double value = 0.0; if (!(line_stream >> value) || !std::isfinite(value)) { - fail(state, "expected a finite number for " + expected); + read_occup_fail(state, "expected a finite number for " + expected); } return value; } -void parse_atom_header(ReaderState& state, - const UnitCell& ucell, - const int expected_type, - const int expected_iat, - const int expected_l, - const int nspin) +static void read_occup_atom_header(OccMatReadState& state, + const UnitCell& ucell, + const int expected_type, + const int expected_iat, + const int expected_l, + const int nspin) { - const std::string line = require_data_line(state, "an atom header"); + const std::string line = read_occup_require_data_line(state, "an atom header"); std::istringstream line_stream(line); std::string element; std::string atom_keyword; @@ -106,36 +108,37 @@ void parse_atom_header(ReaderState& state, if (!(line_stream >> element >> atom_keyword >> atom_index >> l_keyword >> l >> mag_keyword)) { - fail(state, "expected compact header ' Atom L mag ...'"); + read_occup_fail(state, "expected compact header ' Atom L mag ...'"); } if (atom_keyword != "Atom" || l_keyword != "L" || mag_keyword != "mag") { - fail(state, "only the compact 'Atom L mag ...' header is supported"); + read_occup_fail(state, "only the compact 'Atom L mag ...' header is supported"); } if (element != ucell.atoms[expected_type].label) { - fail(state, "element label '" + element + "' does not match the unit cell"); + read_occup_fail(state, "element label '" + element + "' does not match the unit cell"); } if (atom_index != expected_iat + 1) { - fail(state, "atom index does not match the unit-cell order"); + read_occup_fail(state, "atom index does not match the unit-cell order"); } if (l != expected_l) { - fail(state, "L does not match the configured DFT+U channel"); + read_occup_fail(state, "L does not match the configured DFT+U channel"); } const int magnetization_components = (nspin == 4) ? 3 : 1; for (int component = 0; component < magnetization_components; ++component) { - read_finite_double(state, line_stream, "magnetization"); + read_occup_finite_double(state, line_stream, "magnetization"); } - require_line_end(state, line_stream, "the compact atom header"); + read_occup_require_line_end(state, line_stream, "the compact atom header"); } -void parse_spin_nelec_label(ReaderState& state, const int expected_spin) +static void read_occup_spin_nelec_label(OccMatReadState& state, + const int expected_spin) { - const std::string line = require_data_line(state, "a spin block label"); + const std::string line = read_occup_require_data_line(state, "a spin block label"); std::istringstream line_stream(line); std::string spin_keyword; std::string nelec_keyword; @@ -144,16 +147,18 @@ void parse_spin_nelec_label(ReaderState& state, const int expected_spin) || spin_keyword != "spin" || spin != expected_spin || nelec_keyword != "nelec") { - fail(state, "expected compact label 'spin " - + std::to_string(expected_spin) + " nelec '"); + read_occup_fail(state, "expected compact label 'spin " + + std::to_string(expected_spin) + " nelec '"); } - read_finite_double(state, line_stream, "nelec"); - require_line_end(state, line_stream, "the spin block label"); + read_occup_finite_double(state, line_stream, "nelec"); + read_occup_require_line_end(state, line_stream, "the spin block label"); } -void parse_soc_offdiagonal_label(ReaderState& state, const std::string& component) +static void read_occup_soc_offdiagonal_label(OccMatReadState& state, + const std::string& component) { - const std::string line = require_data_line(state, "the spin 12 " + component + " label"); + const std::string line + = read_occup_require_data_line(state, "the spin 12 " + component + " label"); std::istringstream line_stream(line); std::string spin_keyword; std::string block; @@ -162,39 +167,41 @@ void parse_soc_offdiagonal_label(ReaderState& state, const std::string& componen || spin_keyword != "spin" || block != "12" || actual_component != component) { - fail(state, "expected compact label 'spin 12 " + component + "'"); + read_occup_fail(state, "expected compact label 'spin 12 " + component + "'"); } - require_line_end(state, line_stream, "the spin 12 " + component + " label"); + read_occup_require_line_end(state, line_stream, "the spin 12 " + component + " label"); } -std::vector parse_matrix(ReaderState& state, const int size) +static std::vector read_occup_matrix(OccMatReadState& state, + const int size) { std::vector matrix(size * size, 0.0); for (int row = 0; row < size; ++row) { - const std::string line = require_data_line(state, "matrix row"); + const std::string line = read_occup_require_data_line(state, "matrix row"); std::istringstream line_stream(line); for (int column = 0; column < size; ++column) { const int index = row * size + column; - matrix[index] = read_finite_double(state, line_stream, "matrix element"); + matrix[index] + = read_occup_finite_double(state, line_stream, "matrix element"); } - require_line_end(state, line_stream, "a matrix row"); + read_occup_require_line_end(state, line_stream, "a matrix row"); } return matrix; } -void read_collinear_atom(ReaderState& state, - OccupationMatrix& occ, - const int iat, - const int l) +static void read_occup_collinear_atom(OccMatReadState& state, + OccupationMatrix& occ, + const int iat, + const int l) { const int size = 2 * l + 1; std::vector> spin_matrices(2); for (int spin = 0; spin < 2; ++spin) { - parse_spin_nelec_label(state, spin + 1); - spin_matrices[spin] = parse_matrix(state, size); + read_occup_spin_nelec_label(state, spin + 1); + spin_matrices[spin] = read_occup_matrix(state, size); } for (int spin = 0; spin < 2; ++spin) @@ -210,10 +217,10 @@ void read_collinear_atom(ReaderState& state, } } -void store_soc_pauli(OccupationMatrix& occ, - const int iat, - const int l, - const std::vector>& blocks) +static void read_occup_store_soc_pauli(OccupationMatrix& occ, + const int iat, + const int l, + const std::vector>& blocks) { ModuleBase::matrix& matrix = occ.mat(iat, l, 0); const int block_size = static_cast(blocks[0].size()); @@ -226,10 +233,11 @@ void store_soc_pauli(OccupationMatrix& occ, } } -void store_soc_spin_basis(OccupationMatrix& occ, - const int iat, - const int l, - const std::vector>& blocks) +static void read_occup_store_soc_spin_basis( + OccupationMatrix& occ, + const int iat, + const int l, + const std::vector>& blocks) { const int size = 2 * l + 1; for (int row = 0; row < size; ++row) @@ -245,78 +253,74 @@ void store_soc_spin_basis(OccupationMatrix& occ, } } -void read_soc_atom(ReaderState& state, - OccupationMatrix& occ, - const int iat, - const int l, - const DFTU_BASE::OccmatSocLayout soc_layout) +static void read_occup_soc_atom(OccMatReadState& state, + OccupationMatrix& occ, + const int iat, + const int l, + const OccmatSocLayout soc_layout) { const int size = 2 * l + 1; std::vector> blocks(4); - parse_spin_nelec_label(state, 1); - blocks[0] = parse_matrix(state, size); - parse_soc_offdiagonal_label(state, "re"); - blocks[1] = parse_matrix(state, size); - parse_soc_offdiagonal_label(state, "im"); - blocks[2] = parse_matrix(state, size); - parse_spin_nelec_label(state, 2); - blocks[3] = parse_matrix(state, size); - - if (soc_layout == DFTU_BASE::SOC_LAYOUT_PAULI) + read_occup_spin_nelec_label(state, 1); + blocks[0] = read_occup_matrix(state, size); + read_occup_soc_offdiagonal_label(state, "re"); + blocks[1] = read_occup_matrix(state, size); + read_occup_soc_offdiagonal_label(state, "im"); + blocks[2] = read_occup_matrix(state, size); + read_occup_spin_nelec_label(state, 2); + blocks[3] = read_occup_matrix(state, size); + + if (soc_layout == SOC_LAYOUT_PAULI) { - store_soc_pauli(occ, iat, l, blocks); + read_occup_store_soc_pauli(occ, iat, l, blocks); } else { - store_soc_spin_basis(occ, iat, l, blocks); + read_occup_store_soc_spin_basis(occ, iat, l, blocks); } } -void validate_reader_configuration(const ReaderState& state, - const UnitCell& ucell, - const std::vector& l_channel, - const int nspin, - const int npol, - const DFTU_BASE::OccmatSocLayout soc_layout) +static void read_occup_validate_configuration( + const OccMatReadState& state, + const UnitCell& ucell, + const std::vector& l_channel, + const int nspin, + const int npol, + const OccmatSocLayout soc_layout) { if (static_cast(l_channel.size()) != ucell.ntype) { - fail(state, "l_channel size does not match the number of atom types"); + read_occup_fail(state, "l_channel size does not match the number of atom types"); } for (int type = 0; type < ucell.ntype; ++type) { if (l_channel[type] < -1 || l_channel[type] > ucell.atoms[type].nwl) { - fail(state, "a configured DFT+U channel is outside the unit-cell range"); + read_occup_fail(state, "a configured DFT+U channel is outside the unit-cell range"); } } if (nspin != 1 && nspin != 2 && nspin != 4) { - fail(state, "nspin must be 1, 2 or 4"); + read_occup_fail(state, "nspin must be 1, 2 or 4"); } const int expected_npol = (nspin == 4) ? 2 : 1; if (npol != expected_npol) { - fail(state, "npol is inconsistent with nspin"); + read_occup_fail(state, "npol is inconsistent with nspin"); } - if (soc_layout != DFTU_BASE::SOC_LAYOUT_PAULI - && soc_layout != DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL) + if (soc_layout != SOC_LAYOUT_PAULI + && soc_layout != SOC_LAYOUT_SPIN_BASIS_REAL) { - fail(state, "unknown SOC occupation-matrix storage layout"); + read_occup_fail(state, "unknown SOC occupation-matrix storage layout"); } } -void report_open_error(const std::string& filename) +static void read_occup_report_open_error(const std::string& filename) { - ModuleBase::WARNING_QUIT(read_context, + ModuleBase::WARNING_QUIT(read_occup_context, "Can not open compact occupation-matrix file " + filename); } -} // namespace - -namespace DFTU_BASE -{ - void read_occup_m(const UnitCell& ucell, OccupationMatrix& occ, const std::vector& l_channel, @@ -334,14 +338,14 @@ void read_occup_m(const UnitCell& ucell, std::ifstream input(fn.c_str()); if (!input) { - report_open_error(fn); + read_occup_report_open_error(fn); } GlobalV::ofs_running << " DFT+U: read occupation matrix from " << fn << std::endl; std::cout << " DFT+U: read occupation matrix from " << fn << std::endl; - ReaderState state(input, fn); - validate_reader_configuration(state, ucell, l_channel, nspin, npol, soc_layout); + OccMatReadState state(input, fn); + read_occup_validate_configuration(state, ucell, l_channel, nspin, npol, soc_layout); int atom_count = 0; for (int type = 0; type < ucell.ntype; ++type) @@ -354,14 +358,14 @@ void read_occup_m(const UnitCell& ucell, for (int atom = 0; atom < ucell.atoms[type].na; ++atom) { const int iat = ucell.itia2iat(type, atom); - parse_atom_header(state, ucell, type, iat, l, nspin); + read_occup_atom_header(state, ucell, type, iat, l, nspin); if (nspin == 4) { - read_soc_atom(state, occ, iat, l, soc_layout); + read_occup_soc_atom(state, occ, iat, l, soc_layout); } else { - read_collinear_atom(state, occ, iat, l); + read_occup_collinear_atom(state, occ, iat, l); } ++atom_count; } @@ -369,13 +373,13 @@ void read_occup_m(const UnitCell& ucell, if (atom_count == 0) { - fail(state, "the unit cell has no configured DFT+U atom blocks"); + read_occup_fail(state, "the unit cell has no configured DFT+U atom blocks"); } std::string extra_line; - if (next_data_line(state, extra_line)) + if (read_occup_next_data_line(state, extra_line)) { - fail(state, "unexpected content after the final atom block"); + read_occup_fail(state, "unexpected content after the final atom block"); } } From 13249d04e1ed21da6b0a185411a3a0f7e8d41d97 Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Sat, 17 Oct 2026 21:23:23 +0800 Subject: [PATCH 06/13] Put DFT+U occupation reader entry point first --- .../source_pw/module_pwdft/dftu_base_read.cpp | 150 ++++++++++-------- 1 file changed, 88 insertions(+), 62 deletions(-) diff --git a/source/source_pw/module_pwdft/dftu_base_read.cpp b/source/source_pw/module_pwdft/dftu_base_read.cpp index 1874f6ebe9..b7bed8880b 100644 --- a/source/source_pw/module_pwdft/dftu_base_read.cpp +++ b/source/source_pw/module_pwdft/dftu_base_read.cpp @@ -28,6 +28,94 @@ struct OccMatReadState int line_number; }; +[[noreturn]] static void read_occup_fail(const OccMatReadState& state, + const std::string& message); +static bool read_occup_next_data_line(OccMatReadState& state, std::string& line); +static void read_occup_atom_header(OccMatReadState& state, + const UnitCell& ucell, + int expected_type, + int expected_iat, + int expected_l, + int nspin); +static void read_occup_collinear_atom(OccMatReadState& state, + OccupationMatrix& occ, + int iat, + int l); +static void read_occup_soc_atom(OccMatReadState& state, + OccupationMatrix& occ, + int iat, + int l, + OccmatSocLayout soc_layout); +static void read_occup_validate_configuration(const OccMatReadState& state, + const UnitCell& ucell, + const std::vector& l_channel, + int nspin, + int npol, + OccmatSocLayout soc_layout); +static void read_occup_report_open_error(const std::string& filename); + +void read_occup_m(const UnitCell& ucell, + OccupationMatrix& occ, + const std::vector& l_channel, + const std::string& fn, + int nspin, + int npol, + OccmatSocLayout soc_layout) +{ + ModuleBase::TITLE("DFTU_BASE", "read_occup_m"); + if (GlobalV::MY_RANK != 0) + { + return; + } + + std::ifstream input(fn.c_str()); + if (!input) + { + read_occup_report_open_error(fn); + } + + GlobalV::ofs_running << " DFT+U: read occupation matrix from " << fn << std::endl; + std::cout << " DFT+U: read occupation matrix from " << fn << std::endl; + + OccMatReadState state(input, fn); + read_occup_validate_configuration(state, ucell, l_channel, nspin, npol, soc_layout); + + int atom_count = 0; + for (int type = 0; type < ucell.ntype; ++type) + { + const int l = l_channel[type]; + if (l < 0) + { + continue; + } + for (int atom = 0; atom < ucell.atoms[type].na; ++atom) + { + const int iat = ucell.itia2iat(type, atom); + read_occup_atom_header(state, ucell, type, iat, l, nspin); + if (nspin == 4) + { + read_occup_soc_atom(state, occ, iat, l, soc_layout); + } + else + { + read_occup_collinear_atom(state, occ, iat, l); + } + ++atom_count; + } + } + + if (atom_count == 0) + { + read_occup_fail(state, "the unit cell has no configured DFT+U atom blocks"); + } + + std::string extra_line; + if (read_occup_next_data_line(state, extra_line)) + { + read_occup_fail(state, "unexpected content after the final atom block"); + } +} + [[noreturn]] static void read_occup_fail(const OccMatReadState& state, const std::string& message) { @@ -321,66 +409,4 @@ static void read_occup_report_open_error(const std::string& filename) "Can not open compact occupation-matrix file " + filename); } -void read_occup_m(const UnitCell& ucell, - OccupationMatrix& occ, - const std::vector& l_channel, - const std::string& fn, - int nspin, - int npol, - OccmatSocLayout soc_layout) -{ - ModuleBase::TITLE("DFTU_BASE", "read_occup_m"); - if (GlobalV::MY_RANK != 0) - { - return; - } - - std::ifstream input(fn.c_str()); - if (!input) - { - read_occup_report_open_error(fn); - } - - GlobalV::ofs_running << " DFT+U: read occupation matrix from " << fn << std::endl; - std::cout << " DFT+U: read occupation matrix from " << fn << std::endl; - - OccMatReadState state(input, fn); - read_occup_validate_configuration(state, ucell, l_channel, nspin, npol, soc_layout); - - int atom_count = 0; - for (int type = 0; type < ucell.ntype; ++type) - { - const int l = l_channel[type]; - if (l < 0) - { - continue; - } - for (int atom = 0; atom < ucell.atoms[type].na; ++atom) - { - const int iat = ucell.itia2iat(type, atom); - read_occup_atom_header(state, ucell, type, iat, l, nspin); - if (nspin == 4) - { - read_occup_soc_atom(state, occ, iat, l, soc_layout); - } - else - { - read_occup_collinear_atom(state, occ, iat, l); - } - ++atom_count; - } - } - - if (atom_count == 0) - { - read_occup_fail(state, "the unit cell has no configured DFT+U atom blocks"); - } - - std::string extra_line; - if (read_occup_next_data_line(state, extra_line)) - { - read_occup_fail(state, "unexpected content after the final atom block"); - } -} - } // namespace DFTU_BASE From 335338db4d7d406028330eaedd45af7dd81e8794 Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Sat, 17 Oct 2026 21:24:38 +0800 Subject: [PATCH 07/13] Qualify DFT+U SOC layout constants --- source/source_pw/module_pwdft/dftu_base_read.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/source/source_pw/module_pwdft/dftu_base_read.cpp b/source/source_pw/module_pwdft/dftu_base_read.cpp index b7bed8880b..cf88a5bb8e 100644 --- a/source/source_pw/module_pwdft/dftu_base_read.cpp +++ b/source/source_pw/module_pwdft/dftu_base_read.cpp @@ -358,7 +358,7 @@ static void read_occup_soc_atom(OccMatReadState& state, read_occup_spin_nelec_label(state, 2); blocks[3] = read_occup_matrix(state, size); - if (soc_layout == SOC_LAYOUT_PAULI) + if (soc_layout == DFTU_BASE::SOC_LAYOUT_PAULI) { read_occup_store_soc_pauli(occ, iat, l, blocks); } @@ -396,8 +396,8 @@ static void read_occup_validate_configuration( { read_occup_fail(state, "npol is inconsistent with nspin"); } - if (soc_layout != SOC_LAYOUT_PAULI - && soc_layout != SOC_LAYOUT_SPIN_BASIS_REAL) + if (soc_layout != DFTU_BASE::SOC_LAYOUT_PAULI + && soc_layout != DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL) { read_occup_fail(state, "unknown SOC occupation-matrix storage layout"); } From 8e86f9304d1eddf35fc922252913f782bc283f52 Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Sat, 17 Oct 2026 23:30:56 +0800 Subject: [PATCH 08/13] Simplify compact DFT+U occupation reader --- .../source_pw/module_pwdft/dftu_base_read.cpp | 418 ++++++++---------- 1 file changed, 175 insertions(+), 243 deletions(-) diff --git a/source/source_pw/module_pwdft/dftu_base_read.cpp b/source/source_pw/module_pwdft/dftu_base_read.cpp index cf88a5bb8e..85caaebce8 100644 --- a/source/source_pw/module_pwdft/dftu_base_read.cpp +++ b/source/source_pw/module_pwdft/dftu_base_read.cpp @@ -18,41 +18,39 @@ const char* const read_occup_context = "DFTU_BASE::read_occup_m"; struct OccMatReadState { - OccMatReadState(std::ifstream& input, const std::string& filename) - : input(input), filename(filename), line_number(0) + OccMatReadState(std::ifstream& input, + const std::string& filename, + const UnitCell& ucell, + OccupationMatrix& occ, + int nspin, + OccmatSocLayout soc_layout) + : input(input), + filename(filename), + ucell(ucell), + occ(occ), + nspin(nspin), + soc_layout(soc_layout), + line_number(0) { } std::ifstream& input; const std::string& filename; + const UnitCell& ucell; + OccupationMatrix& occ; + int nspin; + OccmatSocLayout soc_layout; int line_number; }; -[[noreturn]] static void read_occup_fail(const OccMatReadState& state, - const std::string& message); -static bool read_occup_next_data_line(OccMatReadState& state, std::string& line); -static void read_occup_atom_header(OccMatReadState& state, - const UnitCell& ucell, - int expected_type, - int expected_iat, - int expected_l, - int nspin); -static void read_occup_collinear_atom(OccMatReadState& state, - OccupationMatrix& occ, - int iat, - int l); -static void read_occup_soc_atom(OccMatReadState& state, - OccupationMatrix& occ, - int iat, - int l, - OccmatSocLayout soc_layout); -static void read_occup_validate_configuration(const OccMatReadState& state, - const UnitCell& ucell, +[[noreturn]] static void read_occup_format_error(const OccMatReadState& state); +static bool read_occup_data_line(OccMatReadState& state, std::string& line); +static void read_occup_validate_configuration(const UnitCell& ucell, const std::vector& l_channel, int nspin, int npol, OccmatSocLayout soc_layout); -static void read_occup_report_open_error(const std::string& filename); +static void read_occup_atom_block(OccMatReadState& state, int type, int iat, int l); void read_occup_m(const UnitCell& ucell, OccupationMatrix& occ, @@ -68,18 +66,19 @@ void read_occup_m(const UnitCell& ucell, return; } + read_occup_validate_configuration(ucell, l_channel, nspin, npol, soc_layout); + std::ifstream input(fn.c_str()); if (!input) { - read_occup_report_open_error(fn); + ModuleBase::WARNING_QUIT(read_occup_context, + "Can not open compact occupation-matrix file " + fn); } GlobalV::ofs_running << " DFT+U: read occupation matrix from " << fn << std::endl; std::cout << " DFT+U: read occupation matrix from " << fn << std::endl; - OccMatReadState state(input, fn); - read_occup_validate_configuration(state, ucell, l_channel, nspin, npol, soc_layout); - + OccMatReadState state(input, fn, ucell, occ, nspin, soc_layout); int atom_count = 0; for (int type = 0; type < ucell.ntype; ++type) { @@ -90,46 +89,31 @@ void read_occup_m(const UnitCell& ucell, } for (int atom = 0; atom < ucell.atoms[type].na; ++atom) { - const int iat = ucell.itia2iat(type, atom); - read_occup_atom_header(state, ucell, type, iat, l, nspin); - if (nspin == 4) - { - read_occup_soc_atom(state, occ, iat, l, soc_layout); - } - else - { - read_occup_collinear_atom(state, occ, iat, l); - } + read_occup_atom_block(state, type, ucell.itia2iat(type, atom), l); ++atom_count; } } - if (atom_count == 0) - { - read_occup_fail(state, "the unit cell has no configured DFT+U atom blocks"); - } - std::string extra_line; - if (read_occup_next_data_line(state, extra_line)) + if (atom_count == 0 || read_occup_data_line(state, extra_line)) { - read_occup_fail(state, "unexpected content after the final atom block"); + read_occup_format_error(state); } } -[[noreturn]] static void read_occup_fail(const OccMatReadState& state, - const std::string& message) +[[noreturn]] static void read_occup_format_error(const OccMatReadState& state) { - std::ostringstream os; - os << "invalid compact occupation-matrix file " << state.filename; + std::ostringstream message; + message << "Invalid compact occupation-matrix format in " << state.filename; if (state.line_number > 0) { - os << " at line " << state.line_number; + message << " at line " << state.line_number; } - os << ": " << message; - ModuleBase::WARNING_QUIT(read_occup_context, os.str()); + message << "."; + ModuleBase::WARNING_QUIT(read_occup_context, message.str()); } -static bool read_occup_next_data_line(OccMatReadState& state, std::string& line) +static bool read_occup_data_line(OccMatReadState& state, std::string& line) { while (std::getline(state.input, line)) { @@ -144,120 +128,114 @@ static bool read_occup_next_data_line(OccMatReadState& state, std::string& line) return false; } -static std::string read_occup_require_data_line(OccMatReadState& state, - const std::string& expected) +static void read_occup_validate_configuration( + const UnitCell& ucell, + const std::vector& l_channel, + const int nspin, + const int npol, + const OccmatSocLayout soc_layout) { - std::string line; - if (!read_occup_next_data_line(state, line)) + bool valid = static_cast(l_channel.size()) == ucell.ntype; + valid = valid && (nspin == 1 || nspin == 2 || nspin == 4); + valid = valid && npol == ((nspin == 4) ? 2 : 1); + valid = valid && (soc_layout == DFTU_BASE::SOC_LAYOUT_PAULI + || soc_layout == DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); + for (int type = 0; valid && type < ucell.ntype; ++type) { - read_occup_fail(state, "unexpected end of file while reading " + expected); + valid = l_channel[type] >= -1 && l_channel[type] <= ucell.atoms[type].nwl; } - return line; -} - -static void read_occup_require_line_end(OccMatReadState& state, - std::istringstream& line_stream, - const std::string& expected) -{ - std::string extra; - if (line_stream >> extra) + if (!valid) { - read_occup_fail(state, "unexpected token '" + extra + "' after " + expected); + ModuleBase::WARNING_QUIT(read_occup_context, + "Invalid DFT+U occupation reader configuration."); } } -static double read_occup_finite_double(OccMatReadState& state, - std::istringstream& line_stream, - const std::string& expected) +static void read_occup_atom_header(OccMatReadState& state, + const int type, + const int iat, + const int l) { - double value = 0.0; - if (!(line_stream >> value) || !std::isfinite(value)) + std::string line; + if (!read_occup_data_line(state, line)) { - read_occup_fail(state, "expected a finite number for " + expected); + read_occup_format_error(state); } - return value; -} -static void read_occup_atom_header(OccMatReadState& state, - const UnitCell& ucell, - const int expected_type, - const int expected_iat, - const int expected_l, - const int nspin) -{ - const std::string line = read_occup_require_data_line(state, "an atom header"); - std::istringstream line_stream(line); + std::istringstream header(line); std::string element; std::string atom_keyword; std::string l_keyword; std::string mag_keyword; int atom_index = 0; - int l = -1; - if (!(line_stream >> element >> atom_keyword >> atom_index >> l_keyword >> l - >> mag_keyword)) + int file_l = -1; + if (!(header >> element >> atom_keyword >> atom_index >> l_keyword >> file_l + >> mag_keyword) + || element != state.ucell.atoms[type].label || atom_keyword != "Atom" + || atom_index != iat + 1 || l_keyword != "L" || file_l != l + || mag_keyword != "mag") { - read_occup_fail(state, "expected compact header ' Atom L mag ...'"); + read_occup_format_error(state); } - if (atom_keyword != "Atom" || l_keyword != "L" || mag_keyword != "mag") - { - read_occup_fail(state, "only the compact 'Atom L mag ...' header is supported"); - } - if (element != ucell.atoms[expected_type].label) - { - read_occup_fail(state, "element label '" + element + "' does not match the unit cell"); - } - if (atom_index != expected_iat + 1) + + const int component_count = (state.nspin == 4) ? 3 : 1; + for (int component = 0; component < component_count; ++component) { - read_occup_fail(state, "atom index does not match the unit-cell order"); + double magnetization = 0.0; + if (!(header >> magnetization) || !std::isfinite(magnetization)) + { + read_occup_format_error(state); + } } - if (l != expected_l) + std::string extra; + if (header >> extra) { - read_occup_fail(state, "L does not match the configured DFT+U channel"); + read_occup_format_error(state); } +} - const int magnetization_components = (nspin == 4) ? 3 : 1; - for (int component = 0; component < magnetization_components; ++component) +static void read_occup_block_label(OccMatReadState& state, const int block) +{ + std::string line; + if (!read_occup_data_line(state, line)) { - read_occup_finite_double(state, line_stream, "magnetization"); + read_occup_format_error(state); } - read_occup_require_line_end(state, line_stream, "the compact atom header"); -} -static void read_occup_spin_nelec_label(OccMatReadState& state, - const int expected_spin) -{ - const std::string line = read_occup_require_data_line(state, "a spin block label"); - std::istringstream line_stream(line); + std::istringstream label(line); std::string spin_keyword; - std::string nelec_keyword; - int spin = 0; - if (!(line_stream >> spin_keyword >> spin >> nelec_keyword) - || spin_keyword != "spin" || spin != expected_spin - || nelec_keyword != "nelec") + const bool offdiagonal = state.nspin == 4 && (block == 1 || block == 2); + if (offdiagonal) { - read_occup_fail(state, "expected compact label 'spin " - + std::to_string(expected_spin) + " nelec '"); + std::string spin_pair; + std::string component; + const char* expected_component = (block == 1) ? "re" : "im"; + if (!(label >> spin_keyword >> spin_pair >> component) + || spin_keyword != "spin" || spin_pair != "12" + || component != expected_component) + { + read_occup_format_error(state); + } + } + else + { + int spin = 0; + std::string nelec_keyword; + double nelec = 0.0; + const int expected_spin = (block == 0) ? 1 : 2; + if (!(label >> spin_keyword >> spin >> nelec_keyword >> nelec) + || spin_keyword != "spin" || spin != expected_spin + || nelec_keyword != "nelec" || !std::isfinite(nelec)) + { + read_occup_format_error(state); + } } - read_occup_finite_double(state, line_stream, "nelec"); - read_occup_require_line_end(state, line_stream, "the spin block label"); -} -static void read_occup_soc_offdiagonal_label(OccMatReadState& state, - const std::string& component) -{ - const std::string line - = read_occup_require_data_line(state, "the spin 12 " + component + " label"); - std::istringstream line_stream(line); - std::string spin_keyword; - std::string block; - std::string actual_component; - if (!(line_stream >> spin_keyword >> block >> actual_component) - || spin_keyword != "spin" || block != "12" - || actual_component != component) + std::string extra; + if (label >> extra) { - read_occup_fail(state, "expected compact label 'spin 12 " + component + "'"); + read_occup_format_error(state); } - read_occup_require_line_end(state, line_stream, "the spin 12 " + component + " label"); } static std::vector read_occup_matrix(OccMatReadState& state, @@ -266,147 +244,101 @@ static std::vector read_occup_matrix(OccMatReadState& state, std::vector matrix(size * size, 0.0); for (int row = 0; row < size; ++row) { - const std::string line = read_occup_require_data_line(state, "matrix row"); - std::istringstream line_stream(line); + std::string line; + if (!read_occup_data_line(state, line)) + { + read_occup_format_error(state); + } + + std::istringstream row_stream(line); for (int column = 0; column < size; ++column) { - const int index = row * size + column; - matrix[index] - = read_occup_finite_double(state, line_stream, "matrix element"); + double value = 0.0; + if (!(row_stream >> value) || !std::isfinite(value)) + { + read_occup_format_error(state); + } + matrix[row * size + column] = value; + } + + std::string extra; + if (row_stream >> extra) + { + read_occup_format_error(state); } - read_occup_require_line_end(state, line_stream, "a matrix row"); } return matrix; } -static void read_occup_collinear_atom(OccMatReadState& state, - OccupationMatrix& occ, - const int iat, - const int l) +static void read_occup_store_blocks(OccMatReadState& state, + const int iat, + const int l, + const std::vector>& blocks) { const int size = 2 * l + 1; - std::vector> spin_matrices(2); - for (int spin = 0; spin < 2; ++spin) - { - read_occup_spin_nelec_label(state, spin + 1); - spin_matrices[spin] = read_occup_matrix(state, size); - } - - for (int spin = 0; spin < 2; ++spin) + if (state.nspin != 4) { - for (int row = 0; row < size; ++row) + for (int spin = 0; spin < 2; ++spin) { - for (int column = 0; column < size; ++column) + for (int row = 0; row < size; ++row) { - const int index = row * size + column; - occ.set(iat, l, spin, row, column, spin_matrices[spin][index]); + for (int column = 0; column < size; ++column) + { + state.occ.set(iat, + l, + spin, + row, + column, + blocks[spin][row * size + column]); + } } } + return; } -} -static void read_occup_store_soc_pauli(OccupationMatrix& occ, - const int iat, - const int l, - const std::vector>& blocks) -{ - ModuleBase::matrix& matrix = occ.mat(iat, l, 0); - const int block_size = static_cast(blocks[0].size()); - for (int index = 0; index < block_size; ++index) + if (state.soc_layout == DFTU_BASE::SOC_LAYOUT_PAULI) { - matrix.c[0 * block_size + index] = blocks[0][index] + blocks[3][index]; - matrix.c[1 * block_size + index] = 2.0 * blocks[1][index]; - matrix.c[2 * block_size + index] = 2.0 * blocks[2][index]; - matrix.c[3 * block_size + index] = blocks[0][index] - blocks[3][index]; + ModuleBase::matrix& matrix = state.occ.mat(iat, l, 0); + const int block_size = size * size; + for (int index = 0; index < block_size; ++index) + { + matrix.c[0 * block_size + index] = blocks[0][index] + blocks[3][index]; + matrix.c[1 * block_size + index] = 2.0 * blocks[1][index]; + matrix.c[2 * block_size + index] = 2.0 * blocks[2][index]; + matrix.c[3 * block_size + index] = blocks[0][index] - blocks[3][index]; + } + return; } -} -static void read_occup_store_soc_spin_basis( - OccupationMatrix& occ, - const int iat, - const int l, - const std::vector>& blocks) -{ - const int size = 2 * l + 1; for (int row = 0; row < size; ++row) { for (int column = 0; column < size; ++column) { const int index = row * size + column; - occ.set(iat, l, 0, row, column, blocks[0][index]); - occ.set(iat, l, 0, row, size + column, blocks[1][index]); - occ.set(iat, l, 0, size + row, column, blocks[1][index]); - occ.set(iat, l, 0, size + row, size + column, blocks[3][index]); + state.occ.set(iat, l, 0, row, column, blocks[0][index]); + state.occ.set(iat, l, 0, row, size + column, blocks[1][index]); + state.occ.set(iat, l, 0, size + row, column, blocks[1][index]); + state.occ.set(iat, l, 0, size + row, size + column, blocks[3][index]); } } } -static void read_occup_soc_atom(OccMatReadState& state, - OccupationMatrix& occ, - const int iat, - const int l, - const OccmatSocLayout soc_layout) +static void read_occup_atom_block(OccMatReadState& state, + const int type, + const int iat, + const int l) { - const int size = 2 * l + 1; - std::vector> blocks(4); - read_occup_spin_nelec_label(state, 1); - blocks[0] = read_occup_matrix(state, size); - read_occup_soc_offdiagonal_label(state, "re"); - blocks[1] = read_occup_matrix(state, size); - read_occup_soc_offdiagonal_label(state, "im"); - blocks[2] = read_occup_matrix(state, size); - read_occup_spin_nelec_label(state, 2); - blocks[3] = read_occup_matrix(state, size); - - if (soc_layout == DFTU_BASE::SOC_LAYOUT_PAULI) - { - read_occup_store_soc_pauli(occ, iat, l, blocks); - } - else - { - read_occup_store_soc_spin_basis(occ, iat, l, blocks); - } -} + read_occup_atom_header(state, type, iat, l); -static void read_occup_validate_configuration( - const OccMatReadState& state, - const UnitCell& ucell, - const std::vector& l_channel, - const int nspin, - const int npol, - const OccmatSocLayout soc_layout) -{ - if (static_cast(l_channel.size()) != ucell.ntype) - { - read_occup_fail(state, "l_channel size does not match the number of atom types"); - } - for (int type = 0; type < ucell.ntype; ++type) - { - if (l_channel[type] < -1 || l_channel[type] > ucell.atoms[type].nwl) - { - read_occup_fail(state, "a configured DFT+U channel is outside the unit-cell range"); - } - } - if (nspin != 1 && nspin != 2 && nspin != 4) - { - read_occup_fail(state, "nspin must be 1, 2 or 4"); - } - const int expected_npol = (nspin == 4) ? 2 : 1; - if (npol != expected_npol) - { - read_occup_fail(state, "npol is inconsistent with nspin"); - } - if (soc_layout != DFTU_BASE::SOC_LAYOUT_PAULI - && soc_layout != DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL) + const int size = 2 * l + 1; + const int block_count = (state.nspin == 4) ? 4 : 2; + std::vector> blocks(block_count); + for (int block = 0; block < block_count; ++block) { - read_occup_fail(state, "unknown SOC occupation-matrix storage layout"); + read_occup_block_label(state, block); + blocks[block] = read_occup_matrix(state, size); } -} - -static void read_occup_report_open_error(const std::string& filename) -{ - ModuleBase::WARNING_QUIT(read_occup_context, - "Can not open compact occupation-matrix file " + filename); + read_occup_store_blocks(state, iat, l, blocks); } } // namespace DFTU_BASE From b0f923102591de439faf481fedbb9adf21adbeb6 Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Sun, 18 Oct 2026 13:16:42 +0800 Subject: [PATCH 09/13] Use flat occupation-matrix writes in DFT+U reader --- .../source_pw/module_pwdft/dftu_base_read.cpp | 52 +++++++++---------- 1 file changed, 25 insertions(+), 27 deletions(-) diff --git a/source/source_pw/module_pwdft/dftu_base_read.cpp b/source/source_pw/module_pwdft/dftu_base_read.cpp index 85caaebce8..8da35722f7 100644 --- a/source/source_pw/module_pwdft/dftu_base_read.cpp +++ b/source/source_pw/module_pwdft/dftu_base_read.cpp @@ -4,6 +4,7 @@ #include "source_base/global_variable.h" #include "source_cell/unitcell.h" +#include #include #include #include @@ -280,47 +281,44 @@ static void read_occup_store_blocks(OccMatReadState& state, { for (int spin = 0; spin < 2; ++spin) { - for (int row = 0; row < size; ++row) - { - for (int column = 0; column < size; ++column) - { - state.occ.set(iat, - l, - spin, - row, - column, - blocks[spin][row * size + column]); - } - } + state.occ.set_flat(iat, l, spin, blocks[spin]); } return; } + const int block_size = size * size; + std::vector matrix(4 * block_size, 0.0); if (state.soc_layout == DFTU_BASE::SOC_LAYOUT_PAULI) { - ModuleBase::matrix& matrix = state.occ.mat(iat, l, 0); - const int block_size = size * size; for (int index = 0; index < block_size; ++index) { - matrix.c[0 * block_size + index] = blocks[0][index] + blocks[3][index]; - matrix.c[1 * block_size + index] = 2.0 * blocks[1][index]; - matrix.c[2 * block_size + index] = 2.0 * blocks[2][index]; - matrix.c[3 * block_size + index] = blocks[0][index] - blocks[3][index]; + matrix[0 * block_size + index] = blocks[0][index] + blocks[3][index]; + matrix[1 * block_size + index] = 2.0 * blocks[1][index]; + matrix[2 * block_size + index] = 2.0 * blocks[2][index]; + matrix[3 * block_size + index] = blocks[0][index] - blocks[3][index]; } - return; } - - for (int row = 0; row < size; ++row) + else { - for (int column = 0; column < size; ++column) + const int matrix_size = 2 * size; + for (int row = 0; row < size; ++row) { - const int index = row * size + column; - state.occ.set(iat, l, 0, row, column, blocks[0][index]); - state.occ.set(iat, l, 0, row, size + column, blocks[1][index]); - state.occ.set(iat, l, 0, size + row, column, blocks[1][index]); - state.occ.set(iat, l, 0, size + row, size + column, blocks[3][index]); + const int block_offset = row * size; + std::copy_n(blocks[0].begin() + block_offset, + size, + matrix.begin() + row * matrix_size); + std::copy_n(blocks[1].begin() + block_offset, + size, + matrix.begin() + row * matrix_size + size); + std::copy_n(blocks[1].begin() + block_offset, + size, + matrix.begin() + (size + row) * matrix_size); + std::copy_n(blocks[3].begin() + block_offset, + size, + matrix.begin() + (size + row) * matrix_size + size); } } + state.occ.set_flat(iat, l, 0, matrix); } static void read_occup_atom_block(OccMatReadState& state, From 3187635a44df940d014db7c61c3e26e62f541dc6 Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Sun, 18 Oct 2026 13:50:20 +0800 Subject: [PATCH 10/13] Rename DFT+U reader implementation --- source/Makefile.Objects | 2 +- source/source_lcao/module_dftu/unittests/CMakeLists.txt | 4 ++-- source/source_pw/module_pwdft/CMakeLists.txt | 2 +- .../{dftu_base_read.cpp => dftu_base_io_read.cpp} | 0 source/source_pw/module_pwdft/test/CMakeLists.txt | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) rename source/source_pw/module_pwdft/{dftu_base_read.cpp => dftu_base_io_read.cpp} (100%) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 9bd8d147d2..bb746cc6fe 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -952,7 +952,7 @@ OBJS_SRCPW=h_ewald_pw.o\ update_cell_pw.o\ dftu_base.o\ dftu_base_io.o\ - dftu_base_read.o\ + dftu_base_io_read.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 56d98da814..6a50a09ffe 100644 --- a/source/source_lcao/module_dftu/unittests/CMakeLists.txt +++ b/source/source_lcao/module_dftu/unittests/CMakeLists.txt @@ -91,7 +91,7 @@ AddTest( ../../../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_read.cpp + ../../../source_pw/module_pwdft/dftu_base_io_read.cpp ../../../source_pw/module_pwdft/yukawa_screening.cpp ../../../source_estate/occ_matrix.cpp ../../../source_estate/occ_mixer.cpp @@ -137,7 +137,7 @@ AddTest( ../../../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_read.cpp + ../../../source_pw/module_pwdft/dftu_base_io_read.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 0e07989054..404fb05b4f 100644 --- a/source/source_pw/module_pwdft/CMakeLists.txt +++ b/source/source_pw/module_pwdft/CMakeLists.txt @@ -31,7 +31,7 @@ endif() list(APPEND objects dftu_base.cpp dftu_base_io.cpp - dftu_base_read.cpp + dftu_base_io_read.cpp dftu_pw.cpp dftu_pw_tools.cpp yukawa_screening.cpp diff --git a/source/source_pw/module_pwdft/dftu_base_read.cpp b/source/source_pw/module_pwdft/dftu_base_io_read.cpp similarity index 100% rename from source/source_pw/module_pwdft/dftu_base_read.cpp rename to source/source_pw/module_pwdft/dftu_base_io_read.cpp diff --git a/source/source_pw/module_pwdft/test/CMakeLists.txt b/source/source_pw/module_pwdft/test/CMakeLists.txt index 51ebf66857..a6cad027a3 100644 --- a/source/source_pw/module_pwdft/test/CMakeLists.txt +++ b/source/source_pw/module_pwdft/test/CMakeLists.txt @@ -28,7 +28,7 @@ AddTest( SOURCES dftu_base_test.cpp ../dftu_base.cpp ../dftu_base_io.cpp - ../dftu_base_read.cpp + ../dftu_base_io_read.cpp ../dftu_pw_tools.cpp ../yukawa_screening.cpp ../../../source_estate/occ_matrix.cpp From 1d7e3ca5757e3ba3a3afedd31b228da8073c2edc Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Sun, 18 Oct 2026 14:00:28 +0800 Subject: [PATCH 11/13] Allow inline comments in compact DFT+U input --- source/source_pw/module_pwdft/dftu_base_io_read.cpp | 5 +++++ .../source_pw/module_pwdft/test/dftu_base_test.cpp | 13 +++++++++++-- 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/source/source_pw/module_pwdft/dftu_base_io_read.cpp b/source/source_pw/module_pwdft/dftu_base_io_read.cpp index 8da35722f7..fba1938410 100644 --- a/source/source_pw/module_pwdft/dftu_base_io_read.cpp +++ b/source/source_pw/module_pwdft/dftu_base_io_read.cpp @@ -124,6 +124,11 @@ static bool read_occup_data_line(OccMatReadState& state, std::string& line) { continue; } + const std::string::size_type comment = line.find('#', first); + if (comment != std::string::npos) + { + line.erase(comment); + } return true; } return false; 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 3d96a1137d..63cc5d8966 100644 --- a/source/source_pw/module_pwdft/test/dftu_base_test.cpp +++ b/source/source_pw/module_pwdft/test/dftu_base_test.cpp @@ -113,7 +113,8 @@ namespace void write_constant_matrix(std::ostream& os, const int size, const double value, - const int first_row_columns) + const int first_row_columns, + const bool inline_comment = false) { for (int row = 0; row < size; ++row) { @@ -124,6 +125,10 @@ void write_constant_matrix(std::ostream& os, { os << " " << value; } + if (inline_comment) + { + os << " # matrix row"; + } os << "\n"; } } @@ -819,7 +824,11 @@ TEST_F(DFTUBaseTest, ReadOccMatCompactHandlesLongComment) std::ofstream ofs(fn); ASSERT_TRUE(ofs.is_open()); ofs << "# " << std::string(300, 'x') << "\n"; - write_compact_collinear_body(ofs); + ofs << "Fe Atom 1 L 2 mag 0.0 # atom header\n"; + ofs << "spin 1 nelec 0.5 # spin header\n"; + write_constant_matrix(ofs, 5, 0.1, 5, true); + ofs << "spin 2 nelec 0.5 # spin header\n"; + write_constant_matrix(ofs, 5, 0.2, 5, true); } const std::vector l_channel = {2}; From 1b16b3b353c4aa1a84dbc1a092418dbed49462a6 Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Sun, 18 Oct 2026 14:01:37 +0800 Subject: [PATCH 12/13] Make malformed DFT+U row tests focused --- source/source_pw/module_pwdft/test/dftu_base_test.cpp | 4 ++++ 1 file changed, 4 insertions(+) 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 63cc5d8966..b1b53d3c66 100644 --- a/source/source_pw/module_pwdft/test/dftu_base_test.cpp +++ b/source/source_pw/module_pwdft/test/dftu_base_test.cpp @@ -909,6 +909,8 @@ TEST_F(DFTUBaseTest, ReadOccMatRejectsShortMatrixRow) ofs << "Fe Atom 1 L 2 mag 0.0\n"; ofs << "spin 1 nelec 0.5\n"; write_constant_matrix(ofs, 5, 0.1, 4); + ofs << "spin 2 nelec 0.5\n"; + write_constant_matrix(ofs, 5, 0.2, 5); } const std::vector l_channel = {2}; OccupationMatrix occ; @@ -933,6 +935,8 @@ TEST_F(DFTUBaseTest, ReadOccMatRejectsLongMatrixRow) ofs << "Fe Atom 1 L 2 mag 0.0\n"; ofs << "spin 1 nelec 0.5\n"; write_constant_matrix(ofs, 5, 0.1, 6); + ofs << "spin 2 nelec 0.5\n"; + write_constant_matrix(ofs, 5, 0.2, 5); } const std::vector l_channel = {2}; OccupationMatrix occ; From abef2d92fde483e4b96a7ffd712b71cd9addbac7 Mon Sep 17 00:00:00 2001 From: lanshuyue <2000011337@stu.pku.edu.cn> Date: Mon, 19 Oct 2026 02:00:05 +0800 Subject: [PATCH 13/13] Test multi-atom DFT+U occupation reading --- .../module_pwdft/test/dftu_base_test.cpp | 78 +++++++++++++++++++ 1 file changed, 78 insertions(+) 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 b1b53d3c66..56f80a8f04 100644 --- a/source/source_pw/module_pwdft/test/dftu_base_test.cpp +++ b/source/source_pw/module_pwdft/test/dftu_base_test.cpp @@ -142,6 +142,18 @@ void write_compact_collinear_body(std::ostream& os) write_constant_matrix(os, 5, 0.2, 5); } +void write_compact_collinear_atom(std::ostream& os, + const int atom_index, + const double spin_up, + const double spin_down) +{ + os << "Fe Atom " << atom_index << " L 2 mag 0.0\n"; + os << "spin 1 nelec 0.5\n"; + write_constant_matrix(os, 5, spin_up, 5); + os << "spin 2 nelec 0.5\n"; + write_constant_matrix(os, 5, spin_down, 5); +} + } // namespace /// After a true -> false re-initialization the Yukawa object must be @@ -845,6 +857,72 @@ TEST_F(DFTUBaseTest, ReadOccMatCompactHandlesLongComment) EXPECT_DOUBLE_EQ(occ.get(0, 2, 1, 4, 4), 0.2); } +TEST_F(DFTUBaseTest, ReadOccMatCompactMultipleAtomsSkipsUncorrelatedSpecies) +{ + UnitCell multi_cell; + multi_cell.ntype = 2; + multi_cell.nat = 4; + multi_cell.atoms = new Atom[multi_cell.ntype]; + multi_cell.set_atom_flag = true; + + multi_cell.atoms[0].label = "O"; + multi_cell.atoms[0].na = 2; + multi_cell.atoms[0].nwl = 1; + multi_cell.atoms[0].nw = 4; + + Atom& fe = multi_cell.atoms[1]; + fe.label = "Fe"; + fe.na = 2; + fe.nwl = 2; + fe.nw = 9; + fe.iw2l.resize(fe.nw); + fe.iw2n.resize(fe.nw); + fe.iw2m.resize(fe.nw); + int iw = 0; + for (int l = 0; l <= fe.nwl; ++l) + { + for (int m = 0; m < 2 * l + 1; ++m) + { + fe.iw2l[iw] = l; + fe.iw2n[iw] = 0; + fe.iw2m[iw] = m; + ++iw; + } + } + + multi_cell.set_iat2itia(); + multi_cell.itia2iat.create(multi_cell.ntype, 2); + multi_cell.itia2iat(0, 0) = 0; + multi_cell.itia2iat(0, 1) = 1; + multi_cell.itia2iat(1, 0) = 2; + multi_cell.itia2iat(1, 1) = 3; + multi_cell.set_iat2iwt(1); + + const std::string fn = testing::TempDir() + "compact_multi_atom.txt"; + { + std::ofstream ofs(fn); + ASSERT_TRUE(ofs.is_open()); + write_compact_collinear_atom(ofs, 3, 0.1, 0.2); + write_compact_collinear_atom(ofs, 4, 0.3, 0.4); + } + + const std::vector l_channel = {-1, 2}; + OccupationMatrix occ; + occ.init(multi_cell, l_channel, 2, 1); + DFTU_BASE::read_occup_m(multi_cell, + occ, + l_channel, + fn, + 2, + 1, + DFTU_BASE::SOC_LAYOUT_PAULI); + + EXPECT_DOUBLE_EQ(occ.get(2, 2, 0, 0, 0), 0.1); + EXPECT_DOUBLE_EQ(occ.get(2, 2, 1, 4, 4), 0.2); + EXPECT_DOUBLE_EQ(occ.get(3, 2, 0, 0, 0), 0.3); + EXPECT_DOUBLE_EQ(occ.get(3, 2, 1, 4, 4), 0.4); +} + TEST_F(DFTUBaseTest, ReadOccMatCompactSocRealSpinBasis) { const std::string fn = testing::TempDir() + "compact_soc_real.txt";