diff --git a/source/source_pw/module_pwdft/dftu_base_io.cpp b/source/source_pw/module_pwdft/dftu_base_io.cpp index b0cbb9e320..c7b1b1971f 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.cpp +++ b/source/source_pw/module_pwdft/dftu_base_io.cpp @@ -308,316 +308,336 @@ 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) +namespace { - ModuleBase::TITLE("DFTU_BASE", "read_occup_m"); - if (GlobalV::MY_RANK != 0) +template +void read_required(std::ifstream& ifs, T& value) +{ + if (!(ifs >> value)) { - return; + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", + "malformed occupation matrix file"); } +} - std::ifstream ifdftu(fn.c_str(), std::ios::in); +void read_legacy_atom_header(std::ifstream& ifs, int& l) +{ + std::string word; + read_required(ifs, word); + if (word.compare(0, 2, "L=") != 0) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); + } + if (word.size() > 2) + { + l = atoi(word.c_str() + 2); + } + else + { + read_required(ifs, l); + } - if (!ifdftu) + read_required(ifs, word); + int zeta = 0; + if (word.compare(0, 8, "ORBITAL=") == 0) { - if (init_occ_mat > 0) + if (word.size() > 8) { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "Can not find the file dm_onsite_ini.txt. Please check your dm_onsite_ini.txt"); + zeta = atoi(word.c_str() + 8); } 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); - } + read_required(ifs, zeta); } - const std::string open_failed_msg = "Can not open " + fn + " file"; - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", open_failed_msg); } + ifs.ignore(150, '\n'); + if (zeta != 0) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", + "only the first radial channel (ORBITAL=0) is supported"); + } +} - 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; +void read_compact_atom_header(std::ifstream& ifs, int& l) +{ + std::string word; + read_required(ifs, word); + if (word != "L") + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); + } + read_required(ifs, l); + ifs.ignore(150, '\n'); +} - while (true) +// Legacy output writes "Atom=1"; older input may also use "Atom= 1". +bool read_atom_index(std::ifstream& ifs, const std::string& word, + int& iat, bool& legacy_layout) +{ + if (word == "Atom") + { + legacy_layout = false; + read_required(ifs, iat); + return true; + } + if (word.compare(0, 5, "Atom=") != 0) + { + return false; + } + legacy_layout = true; + if (word.size() == 5) + { + read_required(ifs, iat); + } + else { - if (!(ifdftu >> word)) + std::istringstream index(word.substr(5)); + if (!(index >> iat) || !index.eof()) { - break; + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "failed to read atom index"); } + } + return true; +} - // Comment lines ("# ...") carry the provenance header and the - // electronic-step marker in the compact layout; skip the line. +// Find the next legacy "Atom=" or compact "Atom" header. +bool read_atom_header(std::ifstream& ifs, int& iat, int& l, bool& legacy_layout) +{ + std::string word; + while (ifs >> word) + { if (word[0] == '#') { - ifdftu.ignore(150, '\n'); + ifs.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) + if (!read_atom_index(ifs, word, iat, legacy_layout)) { continue; } - const bool legacy_layout = legacy_atom; - - ifdftu >> iat; - iat -= 1; - + --iat; if (legacy_layout) { - ifdftu >> word; + read_legacy_atom_header(ifs, l); + } + else + { + read_compact_atom_header(ifs, l); + } + return true; + } + return false; +} - // 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) +int read_collinear_spin_label(std::ifstream& ifs, const int is, const bool legacy_layout) +{ + std::string word; + int spin = is; + if (legacy_layout) + { + read_required(ifs, word); + if (word != "spin=") + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); + } + read_required(ifs, spin); + --spin; + ifs.ignore(150, '\n'); + } + else + { + read_required(ifs, word); + int spin_idx = 0; + read_required(ifs, spin_idx); + if (word != "spin" || spin_idx != is + 1) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG SPIN INDEX IN Plus_U FILE"); + } + read_required(ifs, word); + if (word != "nelec") + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); + } + double nelec = 0.0; + read_required(ifs, nelec); + } + return spin; +} + +void read_collinear_matrix(std::ifstream& ifs, OccupationMatrix& occ, + const int iat, const int l, const bool legacy_layout) +{ + const int nm = 2 * l + 1; + for (int is = 0; is < 2; ++is) + { + const int spin = read_collinear_spin_label(ifs, is, legacy_layout); + for (int m0 = 0; m0 < nm; ++m0) + { + for (int m1 = 0; m1 < nm; ++m1) { - L = atoi(word + 2); + double value = 0.0; + read_required(ifs, value); + occ.set(iat, l, spin, m0, m1, value); } - else + if (legacy_layout) { - ifdftu >> L; + ifs.ignore(150, '\n'); } - 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) +void read_legacy_soc_matrix(std::ifstream& ifs, OccupationMatrix& occ, + const int iat, const int l, const int npol) +{ + const int nm = 2 * l + 1; + 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) { - if (strlen(word) > 8) - { - zeta = atoi(word + 8); - } - else + for (int ipol1 = 0; ipol1 < npol; ++ipol1) { - ifdftu >> zeta; + const int m1_all = m1 + nm * ipol1; + double value = 0.0; + if (!(ifs >> value)) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "failed to read SOC matrix element"); + } + occ.set(iat, l, 0, m0_all, m1_all, value); } } - else - { - zeta = 0; - } - ifdftu.ignore(150, '\n'); + ifs.ignore(150, '\n'); + } + } +} - if (zeta != 0) - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", - "only the first radial channel (ORBITAL=0) is supported"); - } +void read_compact_soc_matrix(std::ifstream& ifs, OccupationMatrix& occ, + const int iat, const int l, const OccmatSocLayout soc_layout) +{ + const int nm = 2 * l + 1; + std::vector uu(nm * nm); + std::vector re(nm * nm); + std::vector im(nm * nm); + std::vector dd(nm * nm); + char word[20]; + + skip_compact_soc_label(ifs, word); + read_matrix_block(ifs, nm, uu.data()); + skip_compact_soc_label(ifs, word); + read_matrix_block(ifs, nm, re.data()); + skip_compact_soc_label(ifs, word); + read_matrix_block(ifs, nm, im.data()); + skip_compact_soc_label(ifs, word); + read_matrix_block(ifs, nm, dd.data()); + + if (soc_layout == SOC_LAYOUT_PAULI) + { + 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 + } + else // SOC_LAYOUT_SPIN_BASIS_REAL cannot store Im(n_ud). + { + for (int m0 = 0; m0 < nm; ++m0) { - // Compact header, e.g. " Fe Atom 1 L 2 mag 0.03..." - ifdftu >> word; - if (strcmp(word, "L") != 0) + for (int m1 = 0; m1 < nm; ++m1) { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); + 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]); } - 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; +void report_occup_open_error(const int init_occ_mat, const std::string& fn, + const std::string& init_chg) +{ + 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); +} + +} // namespace + +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"); - for (int l = 0; l < NL; l++) + if (GlobalV::MY_RANK != 0) + { + return; + } + + std::ifstream ifdftu(fn.c_str(), std::ios::in); + + if (!ifdftu) + { + report_occup_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; + + int iat = 0; + int l = 0; + bool legacy_layout = false; + while (read_atom_header(ifdftu, iat, l, legacy_layout)) + { + const int type = ucell.iat2it[iat]; + const int nl = ucell.atoms[type].nwl + 1; + for (int channel = 0; channel < nl; ++channel) { - if (l != l_channel[T]) + if (channel != l_channel[type]) { 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); - } - } - } - } + read_collinear_matrix(ifdftu, occ, iat, l, legacy_layout); } - else if (nspin == 4) // SOC + else if (nspin == 4) { 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'); - } - } + read_legacy_soc_matrix(ifdftu, occ, iat, l, npol); } 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]); - } - } - } + read_compact_soc_matrix(ifdftu, occ, iat, l, soc_layout); } } } } - - return; } #ifdef __MPI 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..5ae100d369 100644 --- a/source/source_pw/module_pwdft/test/dftu_base_test.cpp +++ b/source/source_pw/module_pwdft/test/dftu_base_test.cpp @@ -780,6 +780,140 @@ TEST_F(DFTUBaseTest, InitBaseReadsOccMatSocPauliRoundtrip) } } +// The old token layout supports both "L=2" and split "L= 2 ORBITAL= 0". +TEST_F(DFTUBaseTest, ReadOccMatLegacyCollinearHeaders) +{ + const std::vector l_channel = {2}; + const std::string dir = testing::TempDir(); + const std::vector headers = { + "Fe 1 Atom=1 L=2 0 0 0\n", + "Fe 1 Atom= 1 L=2 0 0 0\n", + "Fe 1 Atom=1 L= 2 ORBITAL= 0 0 0 0\n" + }; + for (size_t variant = 0; variant < headers.size(); ++variant) + { + const std::string fn = dir + "legacy_collinear_" + std::to_string(variant) + ".txt"; + std::ofstream ofs(fn); + ASSERT_TRUE(ofs.is_open()); + ofs << "# legacy occupation matrix\n" << headers[variant]; + for (int spin = 0; spin < 2; ++spin) + { + ofs << "spin= " << spin + 1 << "\n"; + for (int m0 = 0; m0 < 5; ++m0) + { + for (int m1 = 0; m1 < 5; ++m1) + { + ofs << " " << 100 * spin + 10 * m0 + m1; + } + ofs << "\n"; + } + } + ofs.close(); + + 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); + for (int spin = 0; spin < 2; ++spin) + { + for (int m0 = 0; m0 < 5; ++m0) + { + for (int m1 = 0; m1 < 5; ++m1) + { + EXPECT_DOUBLE_EQ(occ.get(0, 2, spin, m0, m1), + 100 * spin + 10 * m0 + m1); + } + } + } + } +} + +// The legacy SOC file stores the full real spin-basis matrix in interleaved rows. +TEST_F(DFTUBaseTest, ReadOccMatLegacySoc) +{ + const std::vector l_channel = {2}; + const std::string fn = testing::TempDir() + "legacy_soc.txt"; + std::ofstream ofs(fn); + ASSERT_TRUE(ofs.is_open()); + ofs << "Fe 1 Atom=1 L=2 0 0 0\n"; + for (int m0 = 0; m0 < 5; ++m0) + { + for (int ipol0 = 0; ipol0 < 2; ++ipol0) + { + for (int m1 = 0; m1 < 5; ++m1) + { + for (int ipol1 = 0; ipol1 < 2; ++ipol1) + { + ofs << " " << 1000 * ipol0 + 100 * ipol1 + 10 * m0 + m1; + } + } + ofs << "\n"; + } + } + ofs.close(); + + 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); + for (int m0 = 0; m0 < 5; ++m0) + { + for (int ipol0 = 0; ipol0 < 2; ++ipol0) + { + for (int m1 = 0; m1 < 5; ++m1) + { + for (int ipol1 = 0; ipol1 < 2; ++ipol1) + { + EXPECT_DOUBLE_EQ(occ.get(0, 2, 0, m0 + 5 * ipol0, m1 + 5 * ipol1), + 1000 * ipol0 + 100 * ipol1 + 10 * m0 + m1); + } + } + } + } +} + +// The compact SOC path must fill both real spin-off-diagonal blocks and drop Im(n_ud). +TEST_F(DFTUBaseTest, ReadOccMatCompactSocRealSpinBasis) +{ + const std::vector l_channel = {2}; + 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\n"; + const std::vector labels = { + "spin 1 nelec 1", "spin 12 re", "spin 12 im", "spin 2 nelec 1" + }; + for (int block = 0; block < 4; ++block) + { + ofs << labels[block] << "\n"; + for (int m0 = 0; m0 < 5; ++m0) + { + for (int m1 = 0; m1 < 5; ++m1) + { + ofs << " " << 100 * block + 10 * m0 + m1; + } + ofs << "\n"; + } + } + ofs.close(); + + 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); + for (int m0 = 0; m0 < 5; ++m0) + { + for (int m1 = 0; m1 < 5; ++m1) + { + const int offset = 10 * m0 + m1; + EXPECT_DOUBLE_EQ(occ.get(0, 2, 0, m0, m1), offset); + EXPECT_DOUBLE_EQ(occ.get(0, 2, 0, m0, m1 + 5), 100 + offset); + EXPECT_DOUBLE_EQ(occ.get(0, 2, 0, m0 + 5, m1), 100 + offset); + EXPECT_DOUBLE_EQ(occ.get(0, 2, 0, m0 + 5, m1 + 5), 300 + offset); + } + } +} + /// 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