diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 29d0be3276..bb746cc6fe 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -952,6 +952,7 @@ OBJS_SRCPW=h_ewald_pw.o\ update_cell_pw.o\ dftu_base.o\ dftu_base_io.o\ + dftu_base_io_read.o\ dftu_pw.o\ dftu_pw_tools.o\ yukawa_screening.o\ 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_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.cpp b/source/source_pw/module_pwdft/dftu_base.cpp index 58caa93e78..0bb70b999a 100644 --- a/source/source_pw/module_pwdft/dftu_base.cpp +++ b/source/source_pw/module_pwdft/dftu_base.cpp @@ -169,21 +169,20 @@ 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); + 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 @@ -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,11 +210,11 @@ 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, - 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 b0cbb9e320..a8f680b33e 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 @@ -227,41 +226,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 +272,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. /// @@ -773,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"; @@ -985,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; @@ -1299,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 8b24b39eed..babcc76908 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 @@ -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_io_read.cpp new file mode 100644 index 0000000000..fba1938410 --- /dev/null +++ b/source/source_pw/module_pwdft/dftu_base_io_read.cpp @@ -0,0 +1,347 @@ +#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 +#include +#include +#include + +namespace DFTU_BASE +{ + +const char* const read_occup_context = "DFTU_BASE::read_occup_m"; + +struct OccMatReadState +{ + 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_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_atom_block(OccMatReadState& state, int type, int iat, int l); + +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; + } + + read_occup_validate_configuration(ucell, l_channel, nspin, npol, soc_layout); + + std::ifstream input(fn.c_str()); + if (!input) + { + 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, ucell, occ, nspin, 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) + { + read_occup_atom_block(state, type, ucell.itia2iat(type, atom), l); + ++atom_count; + } + } + + std::string extra_line; + if (atom_count == 0 || read_occup_data_line(state, extra_line)) + { + read_occup_format_error(state); + } +} + +[[noreturn]] static void read_occup_format_error(const OccMatReadState& state) +{ + std::ostringstream message; + message << "Invalid compact occupation-matrix format in " << state.filename; + if (state.line_number > 0) + { + message << " at line " << state.line_number; + } + message << "."; + ModuleBase::WARNING_QUIT(read_occup_context, message.str()); +} + +static bool read_occup_data_line(OccMatReadState& 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; + } + const std::string::size_type comment = line.find('#', first); + if (comment != std::string::npos) + { + line.erase(comment); + } + return true; + } + return false; +} + +static void read_occup_validate_configuration( + const UnitCell& ucell, + const std::vector& l_channel, + const int nspin, + const int npol, + const OccmatSocLayout soc_layout) +{ + 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) + { + valid = l_channel[type] >= -1 && l_channel[type] <= ucell.atoms[type].nwl; + } + if (!valid) + { + ModuleBase::WARNING_QUIT(read_occup_context, + "Invalid DFT+U occupation reader configuration."); + } +} + +static void read_occup_atom_header(OccMatReadState& state, + const int type, + const int iat, + const int l) +{ + std::string line; + if (!read_occup_data_line(state, line)) + { + read_occup_format_error(state); + } + + std::istringstream header(line); + std::string element; + std::string atom_keyword; + std::string l_keyword; + std::string mag_keyword; + int atom_index = 0; + 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_format_error(state); + } + + const int component_count = (state.nspin == 4) ? 3 : 1; + for (int component = 0; component < component_count; ++component) + { + double magnetization = 0.0; + if (!(header >> magnetization) || !std::isfinite(magnetization)) + { + read_occup_format_error(state); + } + } + std::string extra; + if (header >> extra) + { + read_occup_format_error(state); + } +} + +static void read_occup_block_label(OccMatReadState& state, const int block) +{ + std::string line; + if (!read_occup_data_line(state, line)) + { + read_occup_format_error(state); + } + + std::istringstream label(line); + std::string spin_keyword; + const bool offdiagonal = state.nspin == 4 && (block == 1 || block == 2); + if (offdiagonal) + { + 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); + } + } + + std::string extra; + if (label >> extra) + { + read_occup_format_error(state); + } +} + +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) + { + 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) + { + 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); + } + } + return matrix; +} + +static void read_occup_store_blocks(OccMatReadState& state, + const int iat, + const int l, + const std::vector>& blocks) +{ + const int size = 2 * l + 1; + if (state.nspin != 4) + { + for (int spin = 0; spin < 2; ++spin) + { + 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) + { + for (int index = 0; index < block_size; ++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]; + } + } + else + { + const int matrix_size = 2 * size; + for (int row = 0; row < size; ++row) + { + 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, + const int type, + const int iat, + const int l) +{ + read_occup_atom_header(state, type, iat, l); + + 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_block_label(state, block); + blocks[block] = read_occup_matrix(state, size); + } + read_occup_store_blocks(state, iat, l, blocks); +} + +} // 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/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 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..b1b53d3c66 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,43 @@ class DFTUBaseTest : public testing::Test } }; +namespace +{ + +void write_constant_matrix(std::ostream& os, + const int size, + const double value, + const int first_row_columns, + const bool inline_comment = false) +{ + 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; + } + if (inline_comment) + { + os << " # matrix row"; + } + 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 +817,191 @@ 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"; + 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}; + OccupationMatrix occ; + occ.init(ucell, l_channel, 2, 1); + DFTU_BASE::read_occup_m(ucell, + occ, + l_channel, + fn, + 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, + fn, + 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, + fn, + 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); + ofs << "spin 2 nelec 0.5\n"; + write_constant_matrix(ofs, 5, 0.2, 5); + } + 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, + fn, + 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); + ofs << "spin 2 nelec 0.5\n"; + write_constant_matrix(ofs, 5, 0.2, 5); + } + 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, + fn, + 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, + fn, + 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, + fn, + 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