diff --git a/AGENTS.md b/AGENTS.md index 93883bfc676..d251dc181dc 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -58,6 +58,23 @@ rules. Read the complete governance document before making or reviewing changes: followed by `f(..., isk, ...)` over inlining the ternary into the call. 15. Do not use `#pragma once` in header files; use traditional `#ifndef`/`#define`/`#endif` include guards instead. + 16. When reading a file, print a one-line confirmation with its full path to + both stdout and the running log; when creating or overwriting a file, + print the same confirmation to the running log. + 17. Do not create or delete files or directories from inside a unit test + via shell calls (`std::system("mkdir")`, `std::system("rm")`); scratch + files shall be written to `testing::TempDir()` rather than the working + directory. Direct library calls (`std::remove`, + `std::filesystem::remove_all`) on test-written files are acceptable. + 18. Do not use non-reentrant C standard library functions that return + pointers to static internal buffers (`std::localtime`, `std::gmtime`, + `std::asctime`, `std::ctime`, `std::strtok`, etc.) in new or modified + code; use the `_r` suffix variants or equivalent thread-safe + alternatives instead. + 19. When reading structured numerical data from user-supplied files, check + the stream state (`fail()`/`eof()`) after extraction and abort with a + diagnostic message on malformed input rather than silently defaulting + to zero. - Use LF line endings for text files. Only `.bat` and `.cmd` files may use CRLF. - Keep source file additions deterministic: update the relevant `CMakeLists.txt` or explain why the file is generated or included indirectly. diff --git a/docs/advanced/elec_properties/band.md b/docs/advanced/elec_properties/band.md index 19bcca6ff11..2a0f7256a47 100644 --- a/docs/advanced/elec_properties/band.md +++ b/docs/advanced/elec_properties/band.md @@ -8,7 +8,7 @@ out_chg 1 ``` With this input parameter, the converged charge density will be output in the files such as `chgs1.cube`, `chgs2.cube`, etc. -Then, one can use the same `STRU` file, pseudopotential files and atomic orbital files (and the local density matrix file dm_onsite.txt if DFT+U is used) to do a non-self-consistent (NSCF) calculation. In this example, the potential is constructed from the ground-state charge density from the proceeding calculation. Now the INPUT file is like: +Then, one can use the same `STRU` file, pseudopotential files and atomic orbital files (and the local density matrix file occ_mat.txt if DFT+U is used) to do a non-self-consistent (NSCF) calculation. In this example, the potential is constructed from the ground-state charge density from the proceeding calculation. Now the INPUT file is like: ``` INPUT_PARAMETERS diff --git a/docs/advanced/elec_properties/dos.md b/docs/advanced/elec_properties/dos.md index d9f78bc4cc1..3b42d6811fd 100644 --- a/docs/advanced/elec_properties/dos.md +++ b/docs/advanced/elec_properties/dos.md @@ -10,7 +10,7 @@ out_chg 1 ``` this will produce the converged charge density, which is contained in the file chg.cube. -Then, use the same `STRU` file, pseudopotential file and atomic orbital file (and the local density matrix file dm_onsite.txt if DFT+U is used) to do a non-self-consistent calculation. In this example, the potential is constructed from the ground-state charge density from the proceeding calculation. Now the INPUT file is like: +Then, use the same `STRU` file, pseudopotential file and atomic orbital file (and the local density matrix file occ_mat.txt if DFT+U is used) to do a non-self-consistent calculation. In this example, the potential is constructed from the ground-state charge density from the proceeding calculation. Now the INPUT file is like: ``` INPUT_PARAMETERS diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp index eac97871445..9d12dfe07dd 100644 --- a/source/source_esolver/esolver_ks.cpp +++ b/source/source_esolver/esolver_ks.cpp @@ -186,7 +186,12 @@ void ESolver_KS::runner(BaseCell& basecell, const int istep) this->niter = this->maxniter; this->diag_ethr = this->inp_->pw_diag_thr; this->scf_nmax_flag = false; // mohan add 2025-09-21 - for (int iter = 1; iter <= this->maxniter; ++iter) + this->scf_rerun_ = false; // clear any stale rerun request before the loop + this->oscillate_esolver = false; // clear any stale oscillation flag before the loop + int iter = 1; + int rerun_count = 0; + const int max_rerun = 100; // safety net: EXX rerun should never exceed this + while (iter <= this->maxniter) { if(iter == this->maxniter) { @@ -212,6 +217,29 @@ void ESolver_KS::runner(BaseCell& basecell, const int istep) } break; } + + // 7) SCF rerun hook. A subclass (e.g. ESolver_KS_PW via EXX) may set + // scf_rerun_ to request a brand-new SCF run starting at iter=1, + // typically after SCF converged but a post-SCF step (EXX exact + // exchange) needs another full SCF. When set, clear the flag and + // restart the loop with iter=1 so every iter-based guard + // (DFT+U occupation-matrix writers, frequency-gated IO, etc.) sees + // a valid 1-based iter. This replaces the old "iter=0 restart + // signal" trick that wrote 0 into the iter reference. + if (this->scf_rerun_) + { + if (++rerun_count > max_rerun) + { + std::cout << " WARNING: SCF rerun limit (" << max_rerun + << ") reached, stopping." << std::endl; + break; + } + this->scf_rerun_ = false; + iter = 1; + continue; + } + + ++iter; } // end scf iterations // 7) after scf diff --git a/source/source_esolver/esolver_ks.h b/source/source_esolver/esolver_ks.h index 2f7055b5d57..29001c4180a 100644 --- a/source/source_esolver/esolver_ks.h +++ b/source/source_esolver/esolver_ks.h @@ -98,6 +98,15 @@ class ESolver_KS : public ESolver_FP bool oscillate_esolver = false; // whether esolver is oscillated bool scf_nmax_flag = false; // whether scf has reached nmax, mohan add 20250921 + /// Set by a subclass (e.g. ESolver_KS_PW via EXX) to request a fresh SCF + /// rerun from iter=1 after the current electronic iteration. ESolver_KS:: + /// runner watches this flag at the bottom of each loop body and, when set, + /// clears it and restarts the SCF loop with iter=1. This replaces the + /// historical "iter=0 restart signal" trick that wrote 0 into the iter + /// reference, which crashed DFT+U occupation-matrix writers that validate + /// iter>=1. + bool scf_rerun_ = false; + /// General EXX info owned by ESolver, initialized from input parameters. General_Exx_Info general_exx_info_; }; diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 79150af093f..76a9bf3b9fc 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -9,6 +9,7 @@ #include "source_lcao/module_deltaspin/deltaspin_lcao.h" #include "source_lcao/setup_dftu_lcao.h" #include "source_pw/module_pwdft/dftu_base.h" // Plus_U_Base (PW and LCAO share it) +#include "source_pw/module_pwdft/dftu_base_io.h" // append_ion_step_snapshot #include "source_hamilt/hs_matrix_k.h" #include "source_estate/module_charge/chg_symm.h" #include "source_estate/module_charge/chg_dmr.h" @@ -581,7 +582,17 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& const std::vector>& dm_vec = this->dmat.dm->get_dmk_vec(); // 1) calculate the local occupation number matrix and energy correction in DFT+U - finish_dftu_lcao(conv_esolver, this->inp_->dft_plus_u, this->inp_->out_chg[0], this->dftu_.get(), ucell, dm_vec, this->kv, this->p_chgmix->get_mixing_beta(), hamilt_lcao, PARAM.globalv.global_out_dir, this->inp_->nspin, PARAM.globalv.npol, PARAM.globalv.gamma_only_local); + const DFTU_BASE::OccmatOutputCfg occmat_cfg{this->inp_->out_freq_ion, + this->inp_->out_freq_elec, + this->inp_->scf_nmax, + this->inp_->out_occ_mat, + this->inp_->dft_plus_u}; + // Read globalv values once into locals so the repeated DFT+U call sites + // below do not each re-enter the global dependency surface. + const std::string& global_out_dir = PARAM.globalv.global_out_dir; + const int npol = PARAM.globalv.npol; + const bool gamma_only_local = PARAM.globalv.gamma_only_local; + finish_dftu_lcao(conv_esolver, this->inp_->dft_plus_u, this->dftu_.get(), ucell, dm_vec, this->kv, this->p_chgmix->get_mixing_beta(), hamilt_lcao, global_out_dir, this->inp_->nspin, npol, gamma_only_local, istep, iter, occmat_cfg); // mohan add 2025-11: push DFT+U energy from Plus_U instance to ElecState. // Covers both dft_plus_u==1 (new method, energy accumulated by DFTU::contributeHR @@ -602,6 +613,39 @@ void ESolver_KS_LCAO::iter_finish(UnitCell& ucell, const int istep, int& // charge mixing is performed, potential is updated, // HF and kS energies are computed, meta-GGA, Jason and restart ESolver_KS::iter_finish(ucell, istep, iter, conv_esolver); + + // overwrite occ_mat.txt with the latest occupation matrix and the + // actual charge-density residual; ESolver_KS::iter_finish() above has + // computed drho for this electronic step. + DFTU_BASE::write_latest_occmat(*this->dftu_, + ucell, + global_out_dir, + this->inp_->nspin, + npol, + istep, + iter, + this->scf_thr, + this->drho, + occmat_cfg, + DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); + + // append the current electronic-step section to occ_matg{#}.txt + // LCAO computes the occupation matrix in every electronic iteration, so + // the snapshot always carries a real matrix. + DFTU_BASE::append_ion_step_snapshot(*this->dftu_, + ucell, + global_out_dir, + this->inp_->nspin, + npol, + istep, + iter, + conv_esolver, + true, + this->scf_thr, + this->drho, + occmat_cfg, + DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); + const bool precision_switched = this->gint_precision_controller_.update_after_iteration(this->drho, this->scf_thr); this->gint_info_->set_exec_precision(this->gint_precision_controller_.current_precision()); if (precision_switched) diff --git a/source/source_esolver/esolver_ks_lcaopw.cpp b/source/source_esolver/esolver_ks_lcaopw.cpp index 2a9968f1737..5c6ba6e9549 100644 --- a/source/source_esolver/esolver_ks_lcaopw.cpp +++ b/source/source_esolver/esolver_ks_lcaopw.cpp @@ -212,7 +212,7 @@ namespace ModuleESolver { // update exx and redo scf XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); - iter = 0; + this->scf_rerun_ = true; std::cout << " Entering 2nd SCF, where EXX is updated" << std::endl; this->two_level_step++; conv_esolver = false; @@ -238,7 +238,7 @@ namespace ModuleESolver gettimeofday(&t_start, nullptr); this->exx_lip->cal_exx(); - iter = 0; + this->scf_rerun_ = true; this->two_level_step++; timeval t_end; diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index 0996349e1d2..7355d60a8c8 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -19,6 +19,7 @@ #include "source_estate/setup_estate_pw.h" // mohan add 20251005 #include "source_hamilt/module_xc/general_exx_info.h" // for General_Exx_Info type used via general_exx_info_ #include "source_io/module_ctrl/ctrl_output_pw.h" // mohan add 20250927 +#include "source_pw/module_pwdft/dftu_base_io.h" // append_ion_step_snapshot #include "source_lcao/module_deltaspin/spin_constrain.h" #include "source_pw/module_pwdft/deltaspin_pw.h" // mohan add 20250309 #include "source_pw/module_pwdft/setup_dftu_pw.h" // mohan add 20250309 @@ -221,6 +222,11 @@ void ESolver_KS_PW::iter_init(UnitCell& ucell, const int istep, const // update local occupations for DFT+U // should before lambda loop in DeltaSpin + const DFTU_BASE::OccmatOutputCfg init_occmat_cfg{this->inp_->out_freq_ion, + this->inp_->out_freq_elec, + this->inp_->scf_nmax, + this->inp_->out_occ_mat, + this->inp_->dft_plus_u}; DFTU_BASE::iter_init_dftu_pw(iter, istep, *this->dftu_, @@ -228,7 +234,12 @@ void ESolver_KS_PW::iter_init(UnitCell& ucell, const int istep, const this->pelec->wg, ucell, this->p_chgmix, - this->kv.isk.data()); + PARAM.globalv.global_out_dir, + init_occmat_cfg, + this->kv.isk.data(), + this->inp_->nspin, + this->inp_->kpar, + PARAM.globalv.npol); // mohan add 2025-11: push DFT+U energy from Plus_U instance to ElecState if (this->inp_->dft_plus_u) @@ -324,12 +335,79 @@ void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int this->ppcell.cal_effective_D(veff, this->pw_rhod, ucell); } - // Handle EXX-related operations after SCF iteration - exx_helper->iter_finish(this->pelec, &this->chr, this->stp.template get_psi_t(), ucell, *this->inp_, conv_esolver, iter); + // Handle EXX-related operations after SCF iteration. + // EXX may override conv_esolver (true -> false) to request an SCF rerun + // after updating the exact-exchange operator. Capture the SCF-converged + // state before the call so we can detect that case and signal + // ESolver_KS::runner via scf_rerun_ to restart the SCF loop at iter=1. + // This replaces the historical iter=0 restart-signal trick (EXX wrote 0 + // into the iter reference), which crashed DFT+U occupation-matrix writers + // validating iter>=1 and broke the 1-based iteration invariant. + const bool scf_converged = conv_esolver; + exx_helper->iter_finish(this->pelec, + &this->chr, + this->stp.template get_psi_t(), + ucell, + *this->inp_, + conv_esolver, + iter); + if (scf_converged && !conv_esolver) + { + this->scf_rerun_ = true; + } // check if oscillate for delta_spin method pw::check_deltaspin_oscillation(iter, this->drho, this->p_chgmix, *this->inp_); + // overwrite occ_mat.txt with the latest occupation matrix and the + // actual charge-density residual of this electronic step. At the very + // first step (istep 0, iter 1) the PW occupation matrix does not exist + // yet (iter_init_dftu_pw returns before calculating it), so skip it and + // keep the previous behavior where occ_mat.txt first appears at iter 2. + const DFTU_BASE::OccmatOutputCfg occmat_cfg{this->inp_->out_freq_ion, + this->inp_->out_freq_elec, + this->inp_->scf_nmax, + this->inp_->out_occ_mat, + this->inp_->dft_plus_u}; + // Read globalv values once into locals so the repeated DFT+U call sites + // below do not each re-enter the global dependency surface. + const std::string& global_out_dir = PARAM.globalv.global_out_dir; + const int npol = PARAM.globalv.npol; + const bool latest_ready = (iter > 1 || istep > 0); + if (latest_ready) + { + DFTU_BASE::write_latest_occmat(*this->dftu_, + ucell, + global_out_dir, + this->inp_->nspin, + npol, + istep, + iter, + this->scf_thr, + this->drho, + occmat_cfg, + DFTU_BASE::SOC_LAYOUT_PAULI); + } + + // append the current electronic-step section to occ_matg{#}.txt + // At istep 0 / iter 1 the PW occupation matrix does not exist yet, + // unless it was loaded from dm_onsite_ini.txt or a previous + // occ_mat.txt; snapshot an "N/A" placeholder in that case. + const bool occmat_ready = latest_ready || this->dftu_->is_occmat_ready(); + DFTU_BASE::append_ion_step_snapshot(*this->dftu_, + ucell, + global_out_dir, + this->inp_->nspin, + npol, + istep, + iter, + conv_esolver, + occmat_ready, + this->scf_thr, + this->drho, + occmat_cfg, + DFTU_BASE::SOC_LAYOUT_PAULI); + // the output quantities ModuleIO::ctrl_iter_pw(istep, iter, conv_esolver, this->stp, this->kv, this->pw_wfc, *this->inp_, GlobalV::ofs_running); } diff --git a/source/source_estate/elecstate_print.cpp b/source/source_estate/elecstate_print.cpp index 4b978c65f70..457c9483c89 100644 --- a/source/source_estate/elecstate_print.cpp +++ b/source/source_estate/elecstate_print.cpp @@ -190,7 +190,7 @@ void print_etot(const Magnetism& magnet, GlobalV::ofs_running << std::setprecision(6); GlobalV::ofs_running << std::setiosflags(std::ios::right); - GlobalV::ofs_running << " Electron density deviation " << scf_thr << std::endl; + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Electron density deviation", scf_thr); if (inp.basis_type == "pw") { @@ -201,7 +201,12 @@ void print_etot(const Magnetism& magnet, std::vector energies_Ry; std::vector energies_eV; - if ((iter % inp.out_freq_elec == 0) || converged || iter == inp.scf_nmax) + // For DFT+U calculations, print the full energy table at every + // electronic step so E_plusU can be monitored. Charge-density and + // wavefunction output remain gated by out_freq_elec elsewhere. + const bool print_table = (iter % inp.out_freq_elec == 0) || converged + || (iter == inp.scf_nmax) || inp.dft_plus_u; + if (print_table) { int n_order = std::max(0, Occupy::gaussian_type); diff --git a/source/source_estate/occ_mixer.h b/source/source_estate/occ_mixer.h index a1d5f94eec1..c31aca8350c 100644 --- a/source/source_estate/occ_mixer.h +++ b/source/source_estate/occ_mixer.h @@ -52,7 +52,7 @@ class OccMatMixer /** * @brief Seed uom_save from an occupation matrix loaded from file. * - * Used when occ_mat_ctrl != 0 (restart from dm_onsite_ini.txt) so that + * Used when init_occ_mat != 0 (restart from dm_onsite_ini.txt) so that * the first mixing step has a meaningful "previous" matrix. */ void seed_save(const OccupationMatrix& occmat); diff --git a/source/source_estate/test/elecstate_print_test.cpp b/source/source_estate/test/elecstate_print_test.cpp index 6746461374e..c2d5bec1743 100644 --- a/source/source_estate/test/elecstate_print_test.cpp +++ b/source/source_estate/test/elecstate_print_test.cpp @@ -188,7 +188,7 @@ TEST_F(ElecStatePrintTest, PrintEtot) GlobalV::ofs_running.close(); ifs.open("test.dat", std::ios::in); std::string str((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); - EXPECT_THAT(str, testing::HasSubstr("Electron density deviation 0.1")); + EXPECT_THAT(str, testing::HasSubstr("Electron density deviation = 0.1")); EXPECT_THAT(str, testing::HasSubstr("Diago Threshold = 0.1")); EXPECT_THAT(str, testing::HasSubstr("E_KohnSham")); EXPECT_THAT(str, testing::HasSubstr("E_vdwD2")); diff --git a/source/source_io/module_parameter/input_conv.cpp b/source/source_io/module_parameter/input_conv.cpp index b05f2826351..8875ecbf764 100644 --- a/source/source_io/module_parameter/input_conv.cpp +++ b/source/source_io/module_parameter/input_conv.cpp @@ -92,7 +92,7 @@ void Input_Conv::Convert() // iteration (1/3) //---------------------------------------------------------- - // Note: DFT+U static members (u_current, u_target, occ_mat_ctrl, etc.) + // Note: DFT+U static members (u_current, u_target, init_occ_mat, etc.) // are now initialized inside Plus_U_Base::init_base() which is called // from setup_pot.cpp (PW) and lcao_set.cpp (LCAO). Mohan refactor 2025-11. diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index a8b584632e0..8f07b087a3f 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -442,6 +442,8 @@ struct Input_para ///< 0: output only when ion steps are finished int out_freq_td = 0; ///< output interval in RT-TDDFT std::vector out_chg = {0, 3}; ///< output charge density. 0: no; 1: yes + bool out_occ_mat = true; ///< output DFT+U occupation matrices; + ///< frequency controlled by out_freq_ion/out_freq_elec std::vector out_xc_r = {-1, 3}; ///< output xc(r). -1: no; >=0: output the order of xc(r) std::vector out_pot = {0, 8}; ///< output potential int out_wfc_pw = 0; ///< 0: no; 1: txt; 2: dat @@ -659,7 +661,10 @@ struct Input_para bool yukawa_potential = false; ///< default: false double yukawa_lambda = -1.0; ///< default: -1.0, which means we calculate lambda double uramping_eV = -1.0; ///< U-Ramping method (eV) - int occ_mat_ctrl = 0; ///< the mode of occupation matrix control + int init_occ_mat = 0; ///< mode of occupation-matrix initialization: 0: no; + ///< 1: read the initial matrix from file, then update it in SCF; + ///< 2: read the matrix from file and keep it fixed + int omc = 0; ///< legacy alias for init_occ_mat double onsite_radius = 0.0; ///< radius of the sphere for onsite projection (Bohr) std::vector hubbard_u_eV = {}; ///< Hubbard Coulomb interaction parameter U(ev) std::vector l_channel = {}; ///< which correlated orbitals need corrected ; d:2 ,f:3, do not 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 7eff18f20a8..bbecf5cace4 100644 --- a/source/source_io/module_parameter/read_inp_exx_dftu.cpp +++ b/source/source_io/module_parameter/read_inp_exx_dftu.cpp @@ -807,19 +807,33 @@ void ReadInput::item_dftu() this->add_item(item); } { - Input_Item item("omc"); - item.annotation = "the mode of occupation matrix control"; + Input_Item item("init_occ_mat"); + item.annotation = "the mode of occupation-matrix initialization"; item.category = "DFT+U correction"; item.type = "Integer"; - item.description = R"(The parameter controls the form of occupation matrix control used. -* 0: No occupation matrix control is performed, and the onsite density matrix will be calculated from wavefunctions in each SCF step. -* 1: The first SCF step will use an initial density matrix read from a file named dm_onsite_ini.txt, but for later steps, the onsite density matrix will be updated. -* 2: The same onsite density matrix from dm_onsite_ini.txt will be used throughout the entire calculation. + item.description = R"(The parameter controls how the DFT+U occupation matrix is initialized. +* 0: No occupation-matrix initialization is performed, and the onsite density matrix will be calculated from wavefunctions in each SCF step. +* 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 easiest way to create dm_onsite_ini.txt is to run a DFT+U calculation with out_chg=1, look for a file named dm_onsite.txt in the OUT.prefix directory, copy and rename it to dm_onsite_ini.txt. The file dm_onsite_ini.txt should be placed in the directory specified by read_file_dir. 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 -> 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.)"; + item.default_value = "0"; + item.unit = ""; + read_sync_int(input.init_occ_mat); + this->add_item(item); + } + { + // note: the variable must be named `item` because read_sync_int() + // refers to it by name inside its macro expansion + Input_Item item("omc"); + item.annotation = "legacy alias for init_occ_mat"; + item.category = "DFT+U correction"; + item.type = "Integer"; + item.description = "Legacy alias for init_occ_mat with the same meaning. If both omc and " + "init_occ_mat are present, init_occ_mat takes precedence."; item.default_value = "0"; item.unit = ""; - read_sync_int(input.occ_mat_ctrl); + read_sync_int(input.omc); this->add_item(item); } { diff --git a/source/source_io/module_parameter/read_inp_out.cpp b/source/source_io/module_parameter/read_inp_out.cpp index e543ebad497..6fec9b3757f 100644 --- a/source/source_io/module_parameter/read_inp_out.cpp +++ b/source/source_io/module_parameter/read_inp_out.cpp @@ -16,6 +16,7 @@ void ReadInput::item_output() item.category = "Output information"; item.type = "Integer"; item.description = "Controls the output interval in ionic steps. When set to a positive integer, information such as charge density, local potential, electrostatic potential, Hamiltonian matrix, overlap matrix, density matrix, Mulliken population analysis, and structure files (STRU{istep} or STRU{istep}.cif, when out_stru is 1 or 2) is printed every n ionic steps." + "\n\nFor DFT+U calculations with out_occ_mat = 1, a separate occupation-matrix file occ_matg{geom_step}.txt is created for every output ionic step. The converged (or final, if not converged) electronic step is always recorded in that file; when out_freq_elec is set smaller than scf_nmax, additional electronic steps are appended as separate sections, each labelled only with the electronic-step number." "\n\n[NOTE] In RT-TDDFT calculations, this parameter is inactive; output frequency is instead controlled by out_freq_td."; item.default_value = "0"; item.unit = ""; @@ -51,7 +52,8 @@ void ReadInput::item_output() item.annotation = "print information every few electronic steps"; item.category = "Output information"; item.type = "Integer"; - item.description = "Output the charge density (only binary format, controlled by out_chg), wavefunction (controlled by out_wfc_pw) per out_freq_elec electronic iterations. Note that they are always output when converged or reach the maximum iterations scf_nmax."; + item.description = "Output the charge density (only binary format, controlled by out_chg), wavefunction (controlled by out_wfc_pw) per out_freq_elec electronic iterations. Note that they are always output when converged or reach the maximum iterations scf_nmax." + "\n\nIn DFT+U calculations with out_freq_ion > 0 and out_occ_mat = 1, this parameter also selects the electronic iterations recorded as separate sections inside occ_matg{geom_step}.txt (the converged or final electronic step is recorded regardless of this parameter)."; item.default_value = "scf_nmax"; item.unit = ""; item.reset_value = [](const Input_Item& item, Parameter& para) { @@ -117,6 +119,18 @@ In molecular dynamics simulations, the output frequency is controlled by out_fre sync_intvec(input.out_chg, 2, 0); this->add_item(item); } + { + Input_Item item("out_occ_mat"); + item.annotation = "output DFT+U occupation matrices"; + item.category = "Output information"; + item.type = "Boolean"; + item.description = "Whether to output the DFT+U occupation matrices. When set to 1 (default), the latest occupation matrix is written to occ_mat.txt at every electronic step, and a numbered file occ_matg{geom_step}.txt is created for every output ionic step selected by out_freq_ion, in which electronic-step sections are appended according to out_freq_elec." + "\n\nThis parameter only takes effect for DFT+U calculations (dft_plus_u > 0) and is independent of out_chg: disabling the charge density output does not disable the occupation-matrix output."; + item.default_value = "True"; + item.unit = ""; + read_sync_bool(input.out_occ_mat); + this->add_item(item); + } { Input_Item item("out_pot"); item.annotation = "output real space potential (with precision 8)"; diff --git a/source/source_io/module_parameter/read_input.cpp b/source/source_io/module_parameter/read_input.cpp index 6257c6dbad3..2de46279d8b 100644 --- a/source/source_io/module_parameter/read_input.cpp +++ b/source/source_io/module_parameter/read_input.cpp @@ -440,6 +440,34 @@ void ReadInput::read_txt_input(Parameter& param, const std::string& filename) } this->normalize_hs_output_options(param); + this->normalize_occ_mat_options(param); +} + +void ReadInput::normalize_occ_mat_options(Parameter& param) +{ + // init_occ_mat is the current name of the former omc parameter; keep + // accepting omc so that INPUT files written by older versions still work. + const auto item_is_read = [this](const std::string& label) { + const auto item = std::find_if( + this->input_lists.begin(), + this->input_lists.end(), + [&label](const std::pair& entry) { return entry.first == label; }); + return item != this->input_lists.end() && item->second.is_read(); + }; + + const bool init_occ_mat_is_read = item_is_read("init_occ_mat"); + const bool omc_is_read = item_is_read("omc"); + if (init_occ_mat_is_read) + { + if (omc_is_read) + { + ModuleBase::WARNING("ReadInput", "both init_occ_mat and omc are set; init_occ_mat takes precedence"); + } + } + else + { + param.input.init_occ_mat = param.input.omc; + } } void ReadInput::normalize_hs_output_options(Parameter& param) diff --git a/source/source_io/module_parameter/read_input.h b/source/source_io/module_parameter/read_input.h index 37484a21247..ac180dcd4ea 100644 --- a/source/source_io/module_parameter/read_input.h +++ b/source/source_io/module_parameter/read_input.h @@ -76,6 +76,13 @@ class ReadInput * does not depend on the order of keywords in the INPUT file. */ void normalize_hs_output_options(Parameter& param); + /** + * @brief Resolve init_occ_mat and its legacy alias omc. + * + * This is called after all INPUT values have been read so that the result + * does not depend on the order of keywords in the INPUT file. + */ + void normalize_occ_mat_options(Parameter& param); /** * @brief write INPUT file of txt format * diff --git a/source/source_io/test/read_input_ptest.cpp b/source/source_io/test/read_input_ptest.cpp index d49c3f3ac74..30398cd7796 100644 --- a/source/source_io/test/read_input_ptest.cpp +++ b/source/source_io/test/read_input_ptest.cpp @@ -361,7 +361,8 @@ TEST_F(InputParaTest, ParaRead) EXPECT_FALSE(param.inp.yukawa_potential); EXPECT_DOUBLE_EQ(param.inp.yukawa_lambda, -1.0); EXPECT_EQ(param.inp.onsite_radius, 0.0); - EXPECT_EQ(param.inp.occ_mat_ctrl, 0); + EXPECT_EQ(param.inp.init_occ_mat, 0); + EXPECT_EQ(param.inp.omc, 0); EXPECT_FALSE(param.inp.dft_plus_dmft); EXPECT_FALSE(param.inp.rpa); EXPECT_FALSE(param.inp.rpa_out_vel); @@ -489,6 +490,52 @@ TEST_F(InputParaTest, ParaRead) EXPECT_DOUBLE_EQ(param.inp.rdmft_power_alpha, 0.656); } +TEST_F(InputParaTest, OccMatLegacyAlias) +{ + // init_occ_mat is the current name of the former omc parameter: omc + // alone must still drive the occupation-matrix initialization, and + // when both keywords are present the new name takes precedence. + const char* input_cases[3] + = {"INPUT_PARAMETERS\n" + "stru_file ./support/STRU\n" + "omc 2\n", + "INPUT_PARAMETERS\n" + "stru_file ./support/STRU\n" + "init_occ_mat 1\n", + "INPUT_PARAMETERS\n" + "stru_file ./support/STRU\n" + "omc 2\n" + "init_occ_mat 1\n"}; + const int expected[3] = {2, 1, 1}; + + for (int icase = 0; icase < 3; ++icase) + { + if (GlobalV::MY_RANK == 0) + { + std::ofstream aliasfile("./INPUT_occ_mat_alias"); + aliasfile << input_cases[icase]; + aliasfile.close(); + } + MPI_Barrier(MPI_COMM_WORLD); + + ModuleIO::ReadInput readinput(GlobalV::MY_RANK); + Parameter param; + readinput.read_parameters(param, "./INPUT_occ_mat_alias"); + EXPECT_EQ(param.inp.init_occ_mat, expected[icase]); + if (icase == 0) + { + EXPECT_EQ(param.inp.omc, 2); + } + + MPI_Barrier(MPI_COMM_WORLD); + if (GlobalV::MY_RANK == 0) + { + std::remove("./INPUT_occ_mat_alias"); + } + MPI_Barrier(MPI_COMM_WORLD); + } +} + TEST_F(InputParaTest, GgaGradAcceptedRange) { ModuleIO::ReadInput readinput(0); diff --git a/source/source_lcao/lcao_set.cpp b/source/source_lcao/lcao_set.cpp index a512a89bf7d..a65df67855b 100644 --- a/source/source_lcao/lcao_set.cpp +++ b/source/source_lcao/lcao_set.cpp @@ -105,8 +105,9 @@ void LCAO_domain::set_pot( inp.device, PARAM.globalv.hubbard_u, PARAM.globalv.uramping, - inp.occ_mat_ctrl, - inp.mixing_dftu); + inp.init_occ_mat, + inp.mixing_dftu, + DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); } //! 4) init exact exchange calculations diff --git a/source/source_lcao/module_dftu/dftu_nao_op.cpp b/source/source_lcao/module_dftu/dftu_nao_op.cpp index b670df11e87..2dc5c9cb6b5 100644 --- a/source/source_lcao/module_dftu/dftu_nao_op.cpp +++ b/source/source_lcao/module_dftu/dftu_nao_op.cpp @@ -74,7 +74,7 @@ hamilt::DFTU_onsite>::DFTU_onsite(HS_Matrix_K* * * For nspin=1: occ is scaled by 0.5 (since only one spin channel computed) * - Subsequent iterations: occ_mat is computed fresh each iteration from updated DMR * - * Case 2: Occ_mat IS ready (is_occmat_ready, i.e., read from dm_onsite.txt file) + * Case 2: Occ_mat IS ready (is_occmat_ready, i.e., read from occ_mat.txt file) * - First electronic iteration: uses pre-read occ_mat directly without DMR calculation * * Skips DMR-based occ calculation entirely * * Reads occ_mat from stored data via get_occ_mat() diff --git a/source/source_lcao/module_dftu/unittests/CMakeLists.txt b/source/source_lcao/module_dftu/unittests/CMakeLists.txt index 1e436770d62..fb7d40cc6c4 100644 --- a/source/source_lcao/module_dftu/unittests/CMakeLists.txt +++ b/source/source_lcao/module_dftu/unittests/CMakeLists.txt @@ -121,6 +121,36 @@ AddTest( ../../../source_hamilt/operator.cpp KEEP_FEATURE_DEFINITIONS __MPI ) + +AddTest( + TARGET MODULE_DFTU_base_io + LIBS parameter psi base device container symmetry + SOURCES test_dftu_base_io.cpp ../dftu_nao_op.cpp ../dftu_nao_adj.cpp ../dftu_nao_pots.cpp ../dftu_nao_fs_r.cpp ../dftu_nao_fs_reduce.cpp ../dftu_nao_for_r.cpp ../dftu_nao_str_r.cpp ../../../source_estate/module_dm/density_matrix.cpp + ../../../source_estate/module_dm/dmr_gamma.cpp + ../../../source_estate/module_dm/dmr_init.cpp + ../../../source_estate/module_dm/dm_setter.cpp + ../../../source_estate/module_dm/dm_getter.cpp + ../../../source_estate/module_dm/dm_tools.cpp + ../../../source_estate/module_dm/dmr_k.cpp + ../../../source_estate/module_dm/dmr_td.cpp + ../../../source_estate/module_dm/dmr_full.cpp + ../../../source_pw/module_pwdft/dftu_base.cpp + ../../../source_pw/module_pwdft/dftu_base_io.cpp + ../../../source_pw/module_pwdft/yukawa_screening.cpp + ../../../source_estate/occ_matrix.cpp + ../../../source_estate/occ_mixer.cpp + ../../../source_hamilt/module_hcontainer/func_folding.cpp + ../../../source_hamilt/module_hcontainer/base_matrix.cpp + ../../../source_hamilt/module_hcontainer/hcontainer.cpp + ../../../source_hamilt/module_hcontainer/atom_pair.cpp + ../../../source_basis/module_ao/parallel_orbitals.cpp + ../../../source_basis/module_ao/orb_atomic_lm.cpp + ../../module_operator_lcao/test/tmp_mocks.cpp + ../../../source_hamilt/operator.cpp + # Reuse the proven MODULE_DFTU_op link closure; MPI-built base/parallel + # objects must see __MPI across TUs. + KEEP_FEATURE_DEFINITIONS __MPI +) endif() endif() diff --git a/source/source_lcao/module_dftu/unittests/test_dftu_base_io.cpp b/source/source_lcao/module_dftu/unittests/test_dftu_base_io.cpp new file mode 100644 index 00000000000..60432cf5908 --- /dev/null +++ b/source/source_lcao/module_dftu/unittests/test_dftu_base_io.cpp @@ -0,0 +1,76 @@ +#include "gtest/gtest.h" + +#include "source_pw/module_pwdft/dftu_base_io.h" + +#include + +#ifdef __MPI +#include "mpi.h" +#endif + +// Focused tests for the occupation-matrix filename helpers and gating logic +// implemented in dftu_base_io.cpp. + +TEST(DFTUBaseIoTest, IonStepGateZeroFreqDisablesNumberedFiles) +{ + const DFTU_BASE::OccmatOutputCfg cfg{0, 10, 10, true, 1}; + EXPECT_FALSE(DFTU_BASE::is_ion_step_output_step(0, cfg)); + EXPECT_FALSE(DFTU_BASE::is_ion_step_output_step(1, cfg)); + EXPECT_FALSE(DFTU_BASE::is_ion_step_output_step(2, cfg)); +} + +TEST(DFTUBaseIoTest, IonStepGateDivisibleIonStepsOnly) +{ + const DFTU_BASE::OccmatOutputCfg cfg{2, 10, 10, true, 1}; + EXPECT_TRUE(DFTU_BASE::is_ion_step_output_step(0, cfg)); + EXPECT_FALSE(DFTU_BASE::is_ion_step_output_step(1, cfg)); + EXPECT_TRUE(DFTU_BASE::is_ion_step_output_step(2, cfg)); + EXPECT_FALSE(DFTU_BASE::is_ion_step_output_step(3, cfg)); + EXPECT_TRUE(DFTU_BASE::is_ion_step_output_step(4, cfg)); +} + +TEST(DFTUBaseIoTest, ElecGatePeriodicTrigger) +{ + const DFTU_BASE::OccmatOutputCfg cfg{1, 2, 10, true, 1}; + EXPECT_FALSE(DFTU_BASE::is_elec_snapshot_trigger(1, false, cfg)); + EXPECT_TRUE(DFTU_BASE::is_elec_snapshot_trigger(2, false, cfg)); + EXPECT_FALSE(DFTU_BASE::is_elec_snapshot_trigger(3, false, cfg)); + EXPECT_TRUE(DFTU_BASE::is_elec_snapshot_trigger(4, false, cfg)); +} + +TEST(DFTUBaseIoTest, ElecGateScfNmaxTrigger) +{ + const DFTU_BASE::OccmatOutputCfg cfg{1, 10, 10, true, 1}; + EXPECT_FALSE(DFTU_BASE::is_elec_snapshot_trigger(1, false, cfg)); + EXPECT_FALSE(DFTU_BASE::is_elec_snapshot_trigger(9, false, cfg)); + EXPECT_TRUE(DFTU_BASE::is_elec_snapshot_trigger(10, false, cfg)); +} + +TEST(DFTUBaseIoTest, ElecGateConvergenceTrigger) +{ + const DFTU_BASE::OccmatOutputCfg cfg{1, 10, 10, true, 1}; + EXPECT_TRUE(DFTU_BASE::is_elec_snapshot_trigger(3, true, cfg)); + EXPECT_TRUE(DFTU_BASE::is_elec_snapshot_trigger(7, true, cfg)); +} + +TEST(DFTUBaseIoTest, IonStepFilenameStartsFromOne) +{ + const std::string fn0 = DFTU_BASE::gen_ion_step_occ_mat_filename("OUT.ABACUS/", 0); + EXPECT_EQ(fn0, "OUT.ABACUS/occ_matg1.txt"); + + const std::string fn4 = DFTU_BASE::gen_ion_step_occ_mat_filename("OUT.ABACUS/", 4); + EXPECT_EQ(fn4, "OUT.ABACUS/occ_matg5.txt"); +} + +int main(int argc, char** argv) +{ +#ifdef __MPI + MPI_Init(&argc, &argv); +#endif + testing::InitGoogleTest(&argc, argv); + const int result = RUN_ALL_TESTS(); +#ifdef __MPI + MPI_Finalize(); +#endif + return result; +} diff --git a/source/source_lcao/setup_dftu_lcao.cpp b/source/source_lcao/setup_dftu_lcao.cpp index 20ed51eaa00..77fe5e0f5d7 100644 --- a/source/source_lcao/setup_dftu_lcao.cpp +++ b/source/source_lcao/setup_dftu_lcao.cpp @@ -41,7 +41,6 @@ void init_dftu_lcao(int dft_plus_u, template void finish_dftu_lcao(const bool conv_esolver, int dft_plus_u, - bool out_chg, void* dftu, const UnitCell& ucell, const std::vector>& dm_vec, @@ -51,7 +50,10 @@ void finish_dftu_lcao(const bool conv_esolver, const std::string& global_out_dir, int nspin, int npol, - const bool gamma_only_local) + const bool gamma_only_local, + int istep, + int iter, + const DFTU_BASE::OccmatOutputCfg& occmat_cfg) { if (!dft_plus_u) { @@ -65,7 +67,7 @@ void finish_dftu_lcao(const bool conv_esolver, /// new DFT+U method calculates energy in Hamiltonian if (dft_plus_u == 2) { - if (dftu_ptr->get_occ_mat_ctrl() != 2) + if (dftu_ptr->get_init_occ_mat() != 2) { const Parallel_Orbitals* pv = hamilt_lcao_ptr->getHR()->get_paraV(); if (pv != nullptr && hamilt_lcao_ptr != nullptr) @@ -80,7 +82,8 @@ void finish_dftu_lcao(const bool conv_esolver, DFTU_LCAO::cal_energy_correction(*dftu_ptr, ucell, PARAM.inp.nspin); } } - DFTU_BASE::output(*dftu_ptr, ucell, out_chg, global_out_dir, nspin, npol); + DFTU_BASE::output(*dftu_ptr, ucell, global_out_dir, nspin, npol, + istep, iter, occmat_cfg, DFTU_BASE::SOC_LAYOUT_SPIN_BASIS_REAL); /// use the converged occupation matrix for next MD/Relax SCF calculation if (conv_esolver) @@ -92,7 +95,6 @@ void finish_dftu_lcao(const bool conv_esolver, /// Template instantiation template void finish_dftu_lcao(const bool conv_esolver, int dft_plus_u, - bool out_chg, void* dftu, const UnitCell& ucell, const std::vector>& dm_vec, @@ -102,11 +104,13 @@ template void finish_dftu_lcao(const bool conv_esolver, const std::string& global_out_dir, int nspin, int npol, - const bool gamma_only_local); + const bool gamma_only_local, + int istep, + int iter, + const DFTU_BASE::OccmatOutputCfg& occmat_cfg); template void finish_dftu_lcao>(const bool conv_esolver, int dft_plus_u, - bool out_chg, void* dftu, const UnitCell& ucell, const std::vector>>& dm_vec, @@ -116,6 +120,9 @@ template void finish_dftu_lcao>(const bool conv_esolver, const std::string& global_out_dir, int nspin, int npol, - const bool gamma_only_local); + const bool gamma_only_local, + int istep, + int iter, + const DFTU_BASE::OccmatOutputCfg& occmat_cfg); } // namespace ModuleESolver diff --git a/source/source_lcao/setup_dftu_lcao.h b/source/source_lcao/setup_dftu_lcao.h index efd0b9ecd3d..ddd1a79bcd9 100644 --- a/source/source_lcao/setup_dftu_lcao.h +++ b/source/source_lcao/setup_dftu_lcao.h @@ -8,6 +8,11 @@ class LCAO_Orbitals; +namespace DFTU_BASE +{ +struct OccmatOutputCfg; // defined in source_pw/module_pwdft/dftu_base_io.h +} + namespace ModuleESolver { @@ -43,22 +48,24 @@ void init_dftu_lcao(int dft_plus_u, * * @param conv_esolver Whether ESolver has converged * @param dft_plus_u DFT+U mode (0=disabled, 1=old, 2=new) - * @param out_chg Whether to output dm_onsite.txt * @param dftu DFT+U object * @param ucell Unit cell * @param dm_vec Density matrix vector * @param kv K-vectors * @param mixing_beta Mixing beta parameter * @param hamilt_lcao Hamiltonian LCAO object - * @param global_out_dir Output directory for dm_onsite.txt + * @param global_out_dir Output directory for occ_mat.txt and occ_matg{#}.txt * @param nspin Number of spin channels (1, 2, or 4) * @param npol Number of polarizations * @param gamma_only_local Whether only the Gamma point is used for LCAO + * @param istep Ionic-step index, starting from 0 + * @param iter Electronic-iteration index, starting from 1 + * @param occmat_cfg Configuration for occupation-matrix output; its out_occ_mat + * member is the master switch of the occ_mat*.txt files */ template void finish_dftu_lcao(const bool conv_esolver, int dft_plus_u, - bool out_chg, void* dftu, const UnitCell& ucell, const std::vector>& dm_vec, @@ -68,7 +75,10 @@ void finish_dftu_lcao(const bool conv_esolver, const std::string& global_out_dir, int nspin, int npol, - const bool gamma_only_local); + const bool gamma_only_local, + int istep, + int iter, + const DFTU_BASE::OccmatOutputCfg& occmat_cfg); } // namespace ModuleESolver diff --git a/source/source_pw/module_pwdft/dftu_base.cpp b/source/source_pw/module_pwdft/dftu_base.cpp index 41f6dd993d3..58caa93e78f 100644 --- a/source/source_pw/module_pwdft/dftu_base.cpp +++ b/source/source_pw/module_pwdft/dftu_base.cpp @@ -41,8 +41,9 @@ void Plus_U_Base::init_base(UnitCell& cell, const std::string& device, const std::vector& hubbard_u, const double uramping, - const int occ_mat_ctrl, - const int mixing_dftu) + const int init_occ_mat, + const int mixing_dftu, + const DFTU_BASE::OccmatSocLayout soc_layout) { ModuleBase::TITLE("Plus_U_Base", "init_base"); @@ -52,7 +53,7 @@ void Plus_U_Base::init_base(UnitCell& cell, this->l_channel = l_channel; this->uramping = uramping; - this->occ_mat_ctrl = occ_mat_ctrl; + this->init_occ_mat = init_occ_mat; this->u_target = hubbard_u; this->u_current = hubbard_u; if (uramping > 0.01) @@ -65,6 +66,17 @@ void Plus_U_Base::init_base(UnitCell& cell, this->energy_u = 0.0; + // The occupation-matrix file is read only once for the whole run. + // occmat_.init() below reallocates and clears the storage on every + // ionic step, so preserve the current in-memory matrix first and + // restore it after init() when the file has already been loaded. + OccupationMatrix occmat_kept; + const bool restore_occmat = (init_occ_mat != 0 && this->occmat_file_loaded_); + if (restore_occmat) + { + occmat_kept = this->occmat_; + } + this->occmat_.init(cell, l_channel, nspin, npol); this->uterm_mat_index.resize(cell.nat); @@ -136,32 +148,76 @@ void Plus_U_Base::init_base(UnitCell& cell, this->yukawa_.reset(); } - if (occ_mat_ctrl != 0) + if (init_occ_mat != 0) { - std::stringstream sst; - sst << global_readin_dir << "dm_onsite_ini.txt"; - DFTU_BASE::read_occup_m(cell, this->occmat_, this->l_channel, this->occ_mat_ctrl, - sst.str(), init_chg, nspin, npol); + if (this->occmat_file_loaded_) + { + // The file was read at the first ionic step; keep using the + // in-memory occupation matrix instead of reading the file + // again. For init_occ_mat=2 this matrix never changes; for + // init_occ_mat=1 it carries the result of the previous + // ionic step. + this->occmat_ = occmat_kept; + this->set_occmat_ready(); + this->occmat_.copy_to_save(cell, this->l_channel); + if (this->has_occ_mixer()) + { + // seed the freshly rebuilt mixer with the current matrix + this->occ_mixer().seed_save(this->occmat_); + } + } + 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. + const std::vector candidates = {"dm_onsite_ini.txt", + "occ_mat.txt", + "dm_onsite.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"); + } + DFTU_BASE::read_occup_m(cell, this->occmat_, this->l_channel, this->init_occ_mat, + readin_fn, init_chg, nspin, npol, soc_layout); #ifdef __MPI - DFTU_BASE::local_occup_bcast(cell, this->occmat_, this->l_channel, nspin, npol); + DFTU_BASE::local_occup_bcast(cell, this->occmat_, this->l_channel, nspin, npol); #endif - this->set_occmat_ready(); - this->occmat_.copy_to_save(cell, this->l_channel); - if (this->has_occ_mixer()) - { - // seed the mixing history with the file-loaded occupation matrix - this->occ_mixer().seed_save(this->occmat_); + this->set_occmat_ready(); + this->occmat_.copy_to_save(cell, this->l_channel); + if (this->has_occ_mixer()) + { + // seed the mixing history with the file-loaded occupation matrix + this->occ_mixer().seed_save(this->occmat_); + } + this->occmat_file_loaded_ = true; } } else { if (init_chg == "file") { - std::stringstream sst; - sst << global_readin_dir << "dm_onsite.txt"; - DFTU_BASE::read_occup_m(cell, this->occmat_, this->l_channel, this->occ_mat_ctrl, - sst.str(), init_chg, nspin, npol); + // 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"}; + 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_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 " + "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); #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.h b/source/source_pw/module_pwdft/dftu_base.h index 0b6e970bae5..af590f53305 100644 --- a/source/source_pw/module_pwdft/dftu_base.h +++ b/source/source_pw/module_pwdft/dftu_base.h @@ -4,6 +4,7 @@ #include "source_base/matrix.h" #include "source_estate/occ_matrix.h" #include "source_estate/occ_mixer.h" +#include "source_pw/module_pwdft/dftu_types.h" #include "source_pw/module_pwdft/yukawa_screening.h" #include @@ -47,8 +48,9 @@ class Plus_U_Base const std::string& device, const std::vector& hubbard_u, const double uramping, - const int occ_mat_ctrl, - const int mixing_dftu); + const int init_occ_mat, + const int mixing_dftu, + const DFTU_BASE::OccmatSocLayout soc_layout); void uramping_update(); bool u_converged(); @@ -64,7 +66,7 @@ class Plus_U_Base const std::vector& get_l_channel_vec() const { return l_channel; } double get_uramping() const { return uramping; } - int get_occ_mat_ctrl() const { return occ_mat_ctrl; } + int get_init_occ_mat() const { return init_occ_mat; } UForm get_form() const { return form; } @@ -145,11 +147,18 @@ class Plus_U_Base std::vector l_channel; double uramping = 0.0; - int occ_mat_ctrl = 0; + int init_occ_mat = 0; // --- Occupation matrices --- OccupationMatrix occmat_; + /// Whether the initial occupation matrix has already been read from + /// file during this run. The file is read exactly once; on later + /// ionic steps the in-memory matrix is preserved (it stays fixed for + /// init_occ_mat=2 and evolves from the previous ionic step for + /// init_occ_mat=1). + bool occmat_file_loaded_ = false; + // Occupation-matrix mixer; constructed only when mixing_dftu != 0. // Owns the flat uom/uom_save buffers and the mixing orchestration. std::unique_ptr occ_mixer_; diff --git a/source/source_pw/module_pwdft/dftu_base_io.cpp b/source/source_pw/module_pwdft/dftu_base_io.cpp index 225b0fdb4c0..b0cbb9e3201 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.cpp +++ b/source/source_pw/module_pwdft/dftu_base_io.cpp @@ -5,14 +5,20 @@ #include "source_base/constants.h" #include "source_base/global_function.h" #include "source_base/global_variable.h" +#include "source_base/module_external/lapack_connector.h" #include "source_base/parallel_common.h" #include "source_base/parallel_global.h" #include "source_base/timer.h" +#include "source_main/version.h" +#include #include +#include #include +#include #include #include +#include #include // local helpers for eigenvalue calculation @@ -86,20 +92,231 @@ inline std::vector CalculateEigenvalues(std::vector> return eigenvalues; } +/// @brief Extract the four Pauli-component blocks of one SOC correlated shell. +/// +/// The blocks, each (2l+1)x(2l+1) stored row-major with index m0*m+m1, are +/// the charge channel and the sigma_x/y/z spin channels. The PW storage keeps +/// the four blocks contiguous; the LCAO storage is a real symmetric matrix in +/// spin basis, whose spin-off-diagonal imaginary channel is not stored. +/// +/// @param dftu DFT+U object holding the occupation matrix +/// @param iat atom index +/// @param l angular-momentum channel +/// @param layout storage layout of the nspin == 4 occupation matrix +/// @param blocks [out] four Pauli-component blocks +inline void extract_soc_pauli_blocks(const Plus_U_Base& dftu, + const int iat, + const int l, + const DFTU_BASE::OccmatSocLayout layout, + std::vector>& blocks) +{ + const int m = 2 * l + 1; + const int m2 = m * m; + blocks.assign(4, std::vector(m2, 0.0)); + + const OccupationMatrix& occmat = dftu.occmat(); + if (layout == DFTU_BASE::SOC_LAYOUT_PAULI) + { + const ModuleBase::matrix& occ = occmat.mat(iat, l, 0); + for (int is = 0; is < 4; ++is) + { + for (int k = 0; k < m2; ++k) + { + blocks[is][k] = occ.c[is * m2 + k]; + } + } + } + else // SOC_LAYOUT_SPIN_BASIS_REAL + { + for (int m0 = 0; m0 < m; ++m0) + { + for (int m1 = 0; m1 < m; ++m1) + { + const int k = m0 * m + m1; + const double suu = occmat.get(iat, l, 0, m0, m1); + const double sud = occmat.get(iat, l, 0, m0, m1 + m); + const double sdu = occmat.get(iat, l, 0, m0 + m, m1); + const double sdd = occmat.get(iat, l, 0, m0 + m, m1 + m); + blocks[0][k] = suu + sdd; + blocks[1][k] = sud + sdu; + // The spin-basis storage drops imaginary parts, so the + // sigma_y channel has no contribution. + blocks[2][k] = 0.0; + blocks[3][k] = suu - sdd; + } + } + } +} + +/// @brief Diagonalize the SOC occupation matrix reconstructed in spinor space. +/// +/// The Hermitian spinor matrix is built from the Pauli blocks as +/// n_upup = (b0 + b3)/2, n_dndn = (b0 - b3)/2, +/// n_updown = (b1 + i*b2)/2. +/// Its 2*(2l+1) real eigenvalues are computed with LAPACK zheev. +/// +/// @param blocks four Pauli-component blocks +/// @param m 2l+1 +/// @return eigenvalues sorted ascending by zheev +inline std::vector calculate_spinor_eigenvalues( + const std::vector>& blocks, + const int m) +{ + const int n = 2 * m; + ModuleBase::ComplexMatrix spinor_n(n, n); + for (int m0 = 0; m0 < m; ++m0) + { + for (int m1 = 0; m1 < m; ++m1) + { + const int k = m0 * m + m1; + const double b0 = blocks[0][k]; + const double b1 = blocks[1][k]; + const double b2 = blocks[2][k]; + const double b3 = blocks[3][k]; + + spinor_n(m0, m1) = std::complex(0.5 * (b0 + b3), 0.0); + spinor_n(m0, m1 + m) = std::complex(0.5 * b1, 0.5 * b2); + spinor_n(m0 + m, m1) = std::complex(0.5 * b1, -0.5 * b2); + spinor_n(m0 + m, m1 + m) = std::complex(0.5 * (b0 - b3), 0.0); + } + } + + std::vector eigenvalues(n, 0.0); + std::vector> work(1, std::complex(0.0, 0.0)); + std::vector rwork(std::max(1, 3 * n - 2), 0.0); + int info = 0; + int lwork = -1; + + LapackConnector::zheev('N', 'U', n, spinor_n, n, eigenvalues.data(), + work.data(), lwork, rwork.data(), &info); + if (info != 0) + { + ModuleBase::WARNING_QUIT("calculate_spinor_eigenvalues", + "LAPACK zheev workspace query failed"); + } + lwork = std::max(1, static_cast(work[0].real())); + work.resize(lwork); + LapackConnector::zheev('N', 'U', n, spinor_n, n, eigenvalues.data(), + work.data(), lwork, rwork.data(), &info); + if (info != 0) + { + ModuleBase::WARNING_QUIT("calculate_spinor_eigenvalues", + "LAPACK zheev diagonalization failed"); + } + return eigenvalues; +} + +} // namespace + + +namespace +{ + +// Write the STRU-style provenance header shared by occ_mat.txt and the +// per-ionic-step snapshot files: ABACUS version, local timestamp and the +// 1-based relaxation step this file belongs to. +void write_provenance_header(std::ostream& os, const int istep) +{ + std::time_t now = std::time(nullptr); + std::tm tm_local; + localtime_r(&now, &tm_local); + char time_buf[64]; + std::strftime(time_buf, sizeof(time_buf), "%Y-%m-%d %H:%M:%S", &tm_local); + os << "# ABACUS version: " << VERSION << "\n"; + os << "# Written at " << time_buf << "\n"; + os << "# RELAX STEP " << istep + 1 << "\n"; +} + +// 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 { +bool is_ion_step_output_step(const int istep, const OccmatOutputCfg& cfg) +{ + if (cfg.out_freq_ion <= 0) + { + return false; + } + return (istep % cfg.out_freq_ion == 0); +} + +bool is_elec_snapshot_trigger(const int iter, + const bool conv_esolver, + const OccmatOutputCfg& cfg) +{ + const bool periodic = (iter % cfg.out_freq_elec == 0); + const bool last_step = (iter == cfg.scf_nmax); + return periodic || last_step || conv_esolver; +} + +std::string gen_ion_step_occ_mat_filename(const std::string& out_dir, const int istep) +{ + std::stringstream ss; + ss << out_dir << "occ_matg" << (istep + 1) << ".txt"; + return ss.str(); +} + +std::string find_first_existing_file(const std::string& dir, + const std::vector& candidates) +{ + for (const std::string& name : candidates) + { + const std::string full = dir + name; + std::ifstream probe(full.c_str()); + if (probe.is_open()) + { + return full; + } + } + return std::string(); +} + void read_occup_m(const UnitCell& ucell, OccupationMatrix& occ, const std::vector& l_channel, - const int occ_mat_ctrl, + const int init_occ_mat, const std::string& fn, const std::string& init_chg, int nspin, - int npol) + int npol, + OccmatSocLayout soc_layout) { ModuleBase::TITLE("DFTU_BASE", "read_occup_m"); @@ -112,7 +329,7 @@ void read_occup_m(const UnitCell& ucell, if (!ifdftu) { - if (occ_mat_ctrl > 0) + 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"); } @@ -120,12 +337,18 @@ void read_occup_m(const UnitCell& ucell, { if (init_chg == "file") { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "Can not find the file dm_onsite.txt. Please do scf calculation first"); + 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); } } - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "Can not open dm_onsite.txt file"); + 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); @@ -137,34 +360,71 @@ void read_occup_m(const UnitCell& ucell, int L = 0; int zeta = 0; - ifdftu.rdstate(); - - while (ifdftu.good()) + while (true) { - ifdftu >> word; - if (ifdftu.eof()) + if (!(ifdftu >> word)) { break; } - if (strcmp("Atom=", word) == 0) + // 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 >> iat; - iat -= 1; ifdftu >> word; - if (strcmp("L=", word) != 0) + // 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"); } - ifdftu >> L; + if (strlen(word) > 2) + { + L = atoi(word + 2); + } + else + { + ifdftu >> L; + } ifdftu >> word; - if (strcmp("ORBITAL=", word) != 0) + // The ORBITAL= token is optional in newly written files; when + // absent the radial channel defaults to 0. + if (strncmp("ORBITAL=", word, 8) == 0) { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); + if (strlen(word) > 8) + { + zeta = atoi(word + 8); + } + else + { + ifdftu >> zeta; + } + } + else + { + zeta = 0; } - ifdftu >> zeta; ifdftu.ignore(150, '\n'); if (zeta != 0) @@ -172,59 +432,109 @@ void read_occup_m(const UnitCell& ucell, 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; + T = ucell.iat2it[iat]; + const int NL = ucell.atoms[T].nwl + 1; - for (int l = 0; l < NL; l++) + for (int l = 0; l < NL; l++) + { + if (l != l_channel[T]) { - if (l != l_channel[T]) - { - continue; - } + continue; + } - if (nspin == 1 || nspin == 2) + const int nm = 2 * L + 1; + + if (nspin == 1 || nspin == 2) + { + for (int is = 0; is < 2; is++) { - for (int is = 0; is < 2; is++) + if (legacy_layout) { ifdftu >> word; - if (strcmp("spin=", word) == 0) + if (strcmp("spin=", word) != 0) { - ifdftu >> spin; - spin -= 1; - ifdftu.ignore(150, '\n'); + 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 < 2 * L + 1; m0++) + double value = 0.0; + for (int m0 = 0; m0 < nm; m0++) + { + for (int m1 = 0; m1 < nm; m1++) { - for (int m1 = 0; m1 < 2 * L + 1; m1++) - { - ifdftu >> value; - occ.set(iat, L, spin, m0, m1, value); - } - ifdftu.ignore(150, '\n'); + 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"); } - else + 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 + } + else if (nspin == 4) // SOC + { + if (legacy_layout) { double value = 0.0; - for (int m0 = 0; m0 < 2 * L + 1; m0++) + for (int m0 = 0; m0 < nm; m0++) { for (int ipol0 = 0; ipol0 < npol; ipol0++) { - const int m0_all = m0 + (2 * L + 1) * ipol0; + const int m0_all = m0 + nm * ipol0; - for (int m1 = 0; m1 < 2 * L + 1; m1++) + for (int m1 = 0; m1 < nm; m1++) { for (int ipol1 = 0; ipol1 < npol; ipol1++) { - int m1_all = m1 + (2 * L + 1) * ipol1; + int m1_all = m1 + nm * ipol1; ifdftu >> value; occ.set(iat, L, 0, m0_all, m1_all, value); } @@ -233,18 +543,77 @@ void read_occup_m(const UnitCell& ucell, } } } - } - } - else - { - ModuleBase::WARNING_QUIT("DFTU_BASE::read_occup_m", "WRONG IN READING LOCAL OCCUPATION NUMBER MATRIX FROM Plus_U FILE"); - } + 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); - ifdftu.rdstate(); + skip_compact_soc_label(ifdftu, word); + read_matrix_block(ifdftu, nm, uu.data()); - if (ifdftu.eof() != 0) - { - break; + 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]); + } + } + } + } + } } } @@ -305,81 +674,266 @@ void local_occup_bcast(const UnitCell& ucell, #endif +void prepare_ion_step_file(const std::string& global_out_dir, + const int istep, + const OccmatOutputCfg& cfg) +{ + // Only the root process creates the per-ionic-step file. + if (GlobalV::MY_RANK != 0) + { + return; + } + + // trunc makes a rerun start from a clean file instead of appending + // snapshots left over from a previous calculation. + const std::string ion_step_fn = gen_ion_step_occ_mat_filename(global_out_dir, istep); + std::ofstream ofs_ion_step; + ofs_ion_step.open(ion_step_fn, std::ios::out | std::ios::trunc); + if (!ofs_ion_step) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::prepare_ion_step_file", + "Can't create per-ionic-step occupation-matrix file"); + } + // Provenance header shared with occ_mat.txt. + write_provenance_header(ofs_ion_step, istep); + ofs_ion_step.close(); + return; +} + + void output(const Plus_U_Base& dftu, const UnitCell& ucell, - bool out_chg, const std::string& global_out_dir, int nspin, - int npol) + int npol, + int istep, + int iter, + const OccmatOutputCfg& cfg, + OccmatSocLayout soc_layout) { ModuleBase::TITLE("DFTU_BASE", "output"); - GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl; - GlobalV::ofs_running << " | #DFT+U INFORMATION# |" << std::endl; - GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl; + if (istep < 0) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::output", "istep must be >= 0"); + } + if (iter < 1) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::output", "iter must be >= 1"); + } + if (cfg.out_freq_ion < 0 || cfg.out_freq_elec < 1 || cfg.scf_nmax < 1) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::output", + "invalid occupation-matrix output frequency configuration"); + } + + // The per-type U values and the per-atom occupation matrices are not + // dumped into the running log here: for large cells (thousands of + // atoms) this floods the log at every electronic step. The dedicated + // occupation-matrix files (occ_mat.txt, occ_matg{#}.txt) gated by + // out_occ_mat carry the same information instead. - for (int T = 0; T < ucell.ntype; T++) + // occ_mat.txt is not written here: its section records drho, which + // is only available after the electronic solve. write_latest_occmat() + // overwrites it at the iter_finish stage. + + // At the first electronic step of an output ionic step, (re)create the + // per-ionic-step file and write its header. The LCAO path enters here + // directly; the PW path additionally prepares the file from + // iter_init_dftu_pw(), because it skips output() at istep 0 / iter 1. + // Electronic-step sections are appended later by append_ion_step_snapshot(). + const bool ion_step_output = is_ion_step_output_step(istep, cfg); + if (cfg.out_occ_mat && ion_step_output && iter == 1) { - const int NL = ucell.atoms[T].nwl + 1; + prepare_ion_step_file(global_out_dir, istep, cfg); + } - for (int L = 0; L < NL; L++) - { - const int N = ucell.atoms[T].l_nchi[L]; + return; +} - if (L >= dftu.get_l_channel(T) && dftu.has_l_channel(T)) - { - if (L != dftu.get_l_channel(T)) - { - continue; - } +// Write one electronic-step section shared by occ_mat.txt and the +// per-ionic-step snapshot files: the step marker, the configured +// charge-density convergence threshold, the actual charge-density +// residual and the full occupation matrices. When occmat_ready is false, +// the matrix body is an "N/A" placeholder. +void write_snapshot_section(std::ostream& os, + const Plus_U_Base& dftu, + const UnitCell& ucell, + const int iter, + const double scf_thr, + const double drho, + const bool occmat_ready, + const int nspin, + const int npol, + const OccmatSocLayout soc_layout) +{ + os << "# Electronic step " << iter << "\n"; + os << "# scf_thr " << std::scientific << std::setprecision(8) << scf_thr << "\n"; + os << "# drho " << std::scientific << std::setprecision(8) << drho << "\n"; - if (!dftu.use_yukawa()) - { - GlobalV::ofs_running << " Type=" << T+1 << " L=" << L << " ORBITAL=" << 0 - << " U=" << dftu.get_u_current(T) * ModuleBase::Ry_to_eV << " eV" << std::endl; - } - else - { - double Ueff = (dftu.yukawa().get_U(T, L) - dftu.yukawa().get_J(T, L)) * ModuleBase::Ry_to_eV; - GlobalV::ofs_running << " Type=" << T+1 << " L=" << L << " ORBITAL=" << 0 - << " U=" << dftu.yukawa().get_U(T, L) * ModuleBase::Ry_to_eV << " eV" - << " J=" << dftu.yukawa().get_J(T, L) * ModuleBase::Ry_to_eV << " eV" - << std::endl; - } - } - } + // 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. + if (!occmat_ready) + { + os << "\n N/A\n"; + return; + } + + // Full occupation matrices together with per-atom magnetism. + write_occup_m(dftu, ucell, os, true, nspin, npol, + OCMAT_FMT_READABLE, soc_layout); +} + +void append_ion_step_snapshot(const Plus_U_Base& dftu, + const UnitCell& ucell, + const std::string& global_out_dir, + int nspin, + int npol, + int istep, + int iter, + bool conv_esolver, + bool occmat_ready, + double scf_thr, + double drho, + const OccmatOutputCfg& cfg, + OccmatSocLayout soc_layout) +{ + ModuleBase::TITLE("DFTU_BASE", "append_ion_step_snapshot"); + + // dft_plus_u <= 0 means the occupation matrix was never computed + // (Plus_U_Base::l_channel stays empty); writing here would dereference + // a null l_channel.data() in write_occup_m. Match the documented contract + // that out_occ_mat only takes effect for DFT+U calculations. Check this + // before validating iter/istep/nspin: when DFT+U is off this function is + // a no-op, so an invalid iter (e.g. from an upstream bug) must not crash + // here. Mirrors the guard-first pattern in write_latest_occmat. + if (!cfg.out_occ_mat || cfg.dft_plus_u <= 0) + { + return; + } + + if (nspin != 1 && nspin != 2 && nspin != 4) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::append_ion_step_snapshot", "nspin must be 1, 2 or 4"); + } + if (istep < 0) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::append_ion_step_snapshot", "istep must be >= 0"); + } + if (iter < 1) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::append_ion_step_snapshot", "iter must be >= 1"); + } + // out_freq_ion == 0 is the valid default (no numbered file) and is + // handled by the gates below; only negative values are an error. + if (cfg.out_freq_ion < 0 || cfg.out_freq_elec < 1 || cfg.scf_nmax < 1) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::append_ion_step_snapshot", + "invalid occupation-matrix output frequency configuration"); } - GlobalV::ofs_running << " Local Occupation Matrices for each atom" << std::endl; - write_occup_m(dftu, ucell, GlobalV::ofs_running, true, nspin, npol); + const bool ion_step_output = is_ion_step_output_step(istep, cfg); + const bool elec_trigger = is_elec_snapshot_trigger(iter, conv_esolver, cfg); + if (!ion_step_output || !elec_trigger || GlobalV::MY_RANK != 0) + { + return; + } - // Write dm_onsite.txt - if (out_chg && GlobalV::MY_RANK == 0) + const std::string ion_step_fn = gen_ion_step_occ_mat_filename(global_out_dir, istep); + + // The provenance header written by prepare_ion_step_file() must not be + // counted as a recorded section: scan the file for an existing section + // marker instead of testing the stream position. + bool section_exists = false; { - std::ofstream ofdftu; - ofdftu.open(global_out_dir + "dm_onsite.txt"); - if (!ofdftu) + std::ifstream checker(ion_step_fn.c_str()); + std::string line; + while (std::getline(checker, line)) { - ModuleBase::WARNING_QUIT("DFTU_BASE::output", "Can't create file dm_onsite.txt"); + if (line.rfind("# Electronic step", 0) == 0) + { + section_exists = true; + break; + } } - write_occup_m(dftu, ucell, ofdftu, false, nspin, npol); - ofdftu.close(); } - GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl; - GlobalV::ofs_running << " | # END DFT+U INFO |" << std::endl; - GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>" << std::endl << std::endl; + std::ofstream ofs_ion_step; + ofs_ion_step.open(ion_step_fn, std::ios::app); + if (!ofs_ion_step) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::append_ion_step_snapshot", + "Can't open per-ionic-step occupation-matrix file"); + } + + // Separate consecutive sections by a blank line; the first section + // directly follows the provenance header. + if (section_exists) + { + ofs_ion_step << "\n"; + } + write_snapshot_section(ofs_ion_step, dftu, ucell, iter, scf_thr, drho, + occmat_ready, nspin, npol, soc_layout); + + ofs_ion_step.close(); + + GlobalV::ofs_running << " Append local occupation matrices of electronic step " << iter + << " to file: " << ion_step_fn << std::endl; + return; +} + + +void write_latest_occmat(const Plus_U_Base& dftu, + const UnitCell& ucell, + const std::string& global_out_dir, + int nspin, + int npol, + int istep, + int iter, + double scf_thr, + double drho, + const OccmatOutputCfg& cfg, + OccmatSocLayout soc_layout) +{ + ModuleBase::TITLE("DFTU_BASE", "write_latest_occmat"); + + if (!cfg.out_occ_mat || cfg.dft_plus_u <= 0 || GlobalV::MY_RANK != 0) + { + return; + } + + const std::string latest_fn = global_out_dir + "occ_mat.txt"; + std::ofstream ofdftu; + ofdftu.open(latest_fn, std::ios::out | std::ios::trunc); + if (!ofdftu) + { + ModuleBase::WARNING_QUIT("DFTU_BASE::write_latest_occmat", "Can't create file occ_mat.txt"); + } + + // occ_mat.txt is a single-section snapshot file: the provenance + // header shared with the g files followed by the same section writer. + write_provenance_header(ofdftu, istep); + // Callers of write_latest_occmat() guarantee the matrix exists: the PW + // esolver gates on latest_ready, the LCAO esolver computes it every + // electronic iteration. + write_snapshot_section(ofdftu, dftu, ucell, iter, scf_thr, drho, + true, nspin, npol, soc_layout); + ofdftu.close(); return; } void write_occup_m(const Plus_U_Base& dftu, const UnitCell& ucell, - std::ofstream& ofs, + std::ostream& ofs, bool diag, int nspin, - int npol) + int npol, + OccmatTextFormat fmt, + OccmatSocLayout soc_layout) { ModuleBase::TITLE("DFTU_BASE", "write_occup_m"); @@ -408,104 +962,370 @@ void write_occup_m(const Plus_U_Base& dftu, continue; } - ofs << "\n Atom= " << iat+1; - ofs << " L= " << l; - ofs << " ORBITAL= " << 0 << std::endl; + const std::string& elem_label = ucell.atoms[T].label; + const int index_in_type = I + 1; + const ModuleBase::Vector3& tau = ucell.atoms[T].tau[I]; + const bool readable = (fmt == OCMAT_FMT_READABLE); + + // In the compact collinear/SOC layout the header carries the + // per-atom magnetism, which is only known after the Pauli + // block traces; that header is written in the spin branches. + const bool defer_header = readable + && (nspin == 1 || nspin == 2 || nspin == 4); + if (!defer_header) + { + if (readable) + { + // Readable header, e.g. " Fe Atom 1 L 2". + // The position is not repeated: it is recorded in + // STRU and constant within one ionic (g) step. + ofs << "\n " << elem_label; + ofs << " Atom " << iat + 1; + ofs << " L " << l << std::endl; + } + else + { + // Legacy token layout parsed by read_occup_m(); the + // " " prefix is skipped when scanning + // for "Atom=". + ofs << "\n" << elem_label << " " << index_in_type; + ofs << " Atom=" << iat + 1; + ofs << " L=" << l; + ofs << std::setprecision(8) << std::fixed + << std::setw(14) << tau.x + << std::setw(14) << tau.y + << std::setw(14) << tau.z << std::endl; + } + } if (nspin == 1 || nspin == 2) { - double sum0[2]; - for (int is = 0; is < 2; is++) + const int nm = 2 * l + 1; + double sum0[2] = {0.0, 0.0}; + + if (readable) { + // Compact snapshot layout, e.g. + // " Fe Atom 1 L 2 mag 0.84215429" + // " spin 1 nelec 1.76153075" + // followed by the matrix rows. if (diag) { - std::vector> A(2 * l + 1, std::vector(2 * l + 1)); - for (int m0 = 0; m0 < 2 * l + 1; m0++) + for (int is = 0; is < 2; is++) { - for (int m1 = 0; m1 < 2 * l + 1; m1++) + for (int m0 = 0; m0 < nm; m0++) { - A[m0][m1] = dftu.occmat().get(iat, l, is, m0, m1); + sum0[is] += dftu.occmat().get(iat, l, is, m0, m0); } } - std::vector eigenvalues = CalculateEigenvalues(A, 2 * l + 1); - sum0[is] = 0.0; - ofs << " Eigenvalues for spin=" << is+1 << std::endl; - ofs << std::setprecision(8) << std::fixed; - for (int i = 0; i < 2 * l + 1; i++) - { - ofs << std::setw(12) << eigenvalues[i]; - sum0[is] += eigenvalues[i]; - } - ofs << std::endl; - ofs << " sum is " << std::setw(12) << sum0[is] << std::endl; } - ofs << " spin= " << is+1 << std::endl; - ofs << std::setprecision(8) << std::fixed; - for (int m0 = 0; m0 < 2 * l + 1; m0++) + ofs << "\n " << elem_label; + ofs << " Atom " << iat + 1; + ofs << " L " << l; + if (diag) { - for (int m1 = 0; m1 < 2 * l + 1; m1++) + ofs << " mag " << std::fixed << std::setprecision(8) + << sum0[0] - sum0[1]; + } + ofs << std::endl; + + ofs << std::fixed << std::setprecision(8); + for (int is = 0; is < 2; is++) + { + ofs << " spin " << is + 1; + if (diag) { - ofs << std::setw(12) - << dftu.occmat().get(iat, l, is, m0, m1); + ofs << " nelec " << sum0[is]; } ofs << std::endl; + for (int m0 = 0; m0 < nm; m0++) + { + for (int m1 = 0; m1 < nm; m1++) + { + // One separating space before every + // value; positive values are 10 chars, + // negative ones 11. + ofs << " " << std::setw(10) + << dftu.occmat().get(iat, l, is, m0, m1); + } + ofs << std::endl; + } } } - if (diag) + else { - ofs << std::setw(12) << std::setprecision(8) - << std::fixed << " Magnetism for atom " << iat+1 << ": " << sum0[0] - sum0[1] - << std::endl; + const std::string trace_label = " Trace (electrons) "; + const std::string matrix_label = " spin= "; + for (int is = 0; is < 2; is++) + { + if (diag) + { + std::vector> A(nm, + std::vector(nm)); + for (int m0 = 0; m0 < nm; m0++) + { + for (int m1 = 0; m1 < nm; m1++) + { + A[m0][m1] = dftu.occmat().get(iat, l, is, m0, m1); + } + } + std::vector eigenvalues = CalculateEigenvalues(A, nm); + sum0[is] = 0.0; + ofs << " Eigenvalues for spin=" << is + 1 << std::endl; + ofs << std::setprecision(8) << std::fixed; + for (int i = 0; i < nm; i++) + { + ofs << std::setw(12) << eigenvalues[i]; + sum0[is] += eigenvalues[i]; + } + ofs << std::endl; + ofs << std::fixed << std::setprecision(8); + ofs << trace_label << sum0[is] << std::endl; + } + ofs << matrix_label << is + 1 << std::endl; + // Matrix elements use fixed-point notation with 8 + // fractional digits (~1e-8 absolute precision). + ofs << std::fixed << std::setprecision(8); + for (int m0 = 0; m0 < nm; m0++) + { + for (int m1 = 0; m1 < nm; m1++) + { + // Legacy keeps its historical + // 12-character column alignment. + ofs << " " << std::setw(11) + << dftu.occmat().get(iat, l, is, m0, m1); + } + ofs << std::endl; + } + } + if (diag) + { + ofs << " mag " << std::fixed << std::setprecision(8) + << sum0[0] - sum0[1] << std::endl; + } } } else if (nspin == 4) // SOC { - if (diag) + const int m = 2 * l + 1; + std::vector> blocks; + extract_soc_pauli_blocks(dftu, iat, l, soc_layout, blocks); + + if (readable) { - double sum0[4]; - std::vector> A(2 * l + 1, std::vector(2 * l + 1)); - int index = 0; - for (int is = 0; is < 4; is++) + // Compact snapshot layout aligned with the + // nspin == 2 style, e.g. + // " Fe Atom 1 L 2 mag mx my mz" + // " spin 1 nelec ..." + // followed by the four spin blocks. + double mag_x = 0.0; + double mag_y = 0.0; + double mag_z = 0.0; + for (int mm = 0; mm < m; ++mm) { - for (int m0 = 0; m0 < 2 * l + 1; m0++) + const int k = mm * m + mm; + mag_x += blocks[1][k]; + mag_y += blocks[2][k]; + mag_z += blocks[3][k]; + } + // Spin-resolved electron counts: traces of the + // up/up and down/down blocks. + double n_up = 0.0; + double n_down = 0.0; + for (int mm = 0; mm < m; ++mm) + { + const int k = mm * m + mm; + n_up += 0.5 * (blocks[0][k] + blocks[3][k]); + n_down += 0.5 * (blocks[0][k] - blocks[3][k]); + } + + ofs << "\n " << elem_label; + ofs << " Atom " << iat + 1; + ofs << " L " << l; + if (diag) + { + ofs << " mag " << std::fixed << std::setprecision(8) + << mag_x << " " << mag_y << " " << mag_z; + } + ofs << std::endl; + + ofs << std::fixed << std::setprecision(8); + + // spin 1: spin up/up block (Re). + ofs << " spin 1"; + if (diag) + { + ofs << " nelec " << n_up; + } + ofs << std::endl; + for (int m0 = 0; m0 < m; ++m0) + { + for (int m1 = 0; m1 < m; ++m1) { - for (int m1 = 0; m1 < 2 * l + 1; m1++) - { - A[m0][m1] = dftu.occmat().get(iat, l, 0, m0, m1); - index++; - } + const int k = m0 * m + m1; + const double val = 0.5 * (blocks[0][k] + blocks[3][k]); + ofs << " " << std::setw(10) << val; } - std::vector eigenvalues = CalculateEigenvalues(A, 2 * l + 1); - sum0[is] = 0.0; - ofs << " Eigenvalues for is=" << is << std::endl; - ofs << std::setprecision(8) << std::fixed; - for (int i = 0; i < 2 * l + 1; i++) + ofs << std::endl; + } + + // Off-diagonal Re block. + ofs << " spin 12 re" << std::endl; + for (int m0 = 0; m0 < m; ++m0) + { + for (int m1 = 0; m1 < m; ++m1) { - ofs << std::setw(12) << eigenvalues[i]; - sum0[is] += eigenvalues[i]; + const int k = m0 * m + m1; + const double val = 0.5 * blocks[1][k]; + ofs << " " << std::setw(10) << val; + } + ofs << std::endl; + } + + // Off-diagonal Im block. + ofs << " spin 12 im" << std::endl; + for (int m0 = 0; m0 < m; ++m0) + { + for (int m1 = 0; m1 < m; ++m1) + { + const int k = m0 * m + m1; + const double val = 0.5 * blocks[2][k]; + ofs << " " << std::setw(10) << val; + } + ofs << std::endl; + } + + // spin 2: spin down/down block (Re). + ofs << " spin 2"; + if (diag) + { + ofs << " nelec " << n_down; + } + ofs << std::endl; + for (int m0 = 0; m0 < m; ++m0) + { + for (int m1 = 0; m1 < m; ++m1) + { + const int k = m0 * m + m1; + const double val = 0.5 * (blocks[0][k] - blocks[3][k]); + ofs << " " << std::setw(10) << val; } ofs << std::endl; - ofs << " sum is " << std::setw(12) << sum0[is] << std::endl; } - ofs << std::setw(12) << std::setprecision(8) - << std::fixed << " Magnetism for atom " << iat + 1 << ": " - << sum0[1] << " " << sum0[2] << " " << sum0[3] << std::endl; } - else + else if (diag) { - for (int m0 = 0; m0 < 2 * l + 1; m0++) + // Occupation numbers of the correlated spinor shell. + std::vector eigenvalues = calculate_spinor_eigenvalues(blocks, m); + ofs << " Eigenvalues in spinor space:" << std::endl; + ofs << std::setprecision(8) << std::fixed; + for (int i = 0; i < 2 * m; ++i) { - for (int ipol0 = 0; ipol0 < npol; ipol0++) + ofs << std::setw(12) << eigenvalues[i]; + if ((i + 1) % m == 0) { - const int m0_all = m0 + (2 * l + 1) * ipol0; + ofs << std::endl; + } + } - for (int m1 = 0; m1 < 2 * l + 1; m1++) + // Per-atom charge and spin-moment traces of the shell. + double n_tot = 0.0; + double mag_x = 0.0; + double mag_y = 0.0; + double mag_z = 0.0; + for (int mm = 0; mm < m; ++mm) + { + const int k = mm * m + mm; + n_tot += blocks[0][k]; + mag_x += blocks[1][k]; + mag_y += blocks[2][k]; + mag_z += blocks[3][k]; + } + ofs << std::fixed << std::setprecision(8); + ofs << " Trace (electrons) " << n_tot << std::endl; + ofs << " Magnetism for atom " << iat + 1 << " (mx, my, mz):" + << std::setw(12) << mag_x << std::setw(12) << mag_y + << std::setw(12) << mag_z << std::endl; + + // Matrix elements use fixed-point notation with 8 + // fractional digits (~1e-8 absolute precision). + ofs << std::fixed << std::setprecision(8); + + // Hermitian occupation matrix shown spin block by spin block. + ofs << " Occupation matrix, spin up/up (Re):" << std::endl; + for (int m0 = 0; m0 < m; ++m0) + { + for (int m1 = 0; m1 < m; ++m1) + { + const int k = m0 * m + m1; + const double val = 0.5 * (blocks[0][k] + blocks[3][k]); + ofs << " " << std::setw(10) << val; + } + ofs << std::endl; + } + ofs << " Occupation matrix, spin up/down (Re):" << std::endl; + for (int m0 = 0; m0 < m; ++m0) + { + for (int m1 = 0; m1 < m; ++m1) + { + const int k = m0 * m + m1; + const double val = 0.5 * blocks[1][k]; + ofs << " " << std::setw(10) << val; + } + ofs << std::endl; + } + ofs << " Occupation matrix, spin up/down (Im):" << std::endl; + for (int m0 = 0; m0 < m; ++m0) + { + for (int m1 = 0; m1 < m; ++m1) + { + const int k = m0 * m + m1; + const double val = 0.5 * blocks[2][k]; + ofs << " " << std::setw(10) << val; + } + ofs << std::endl; + } + ofs << " Occupation matrix, spin down/down (Re):" << std::endl; + for (int m0 = 0; m0 < m; ++m0) + { + for (int m1 = 0; m1 < m; ++m1) + { + const int k = m0 * m + m1; + const double val = 0.5 * (blocks[0][k] - blocks[3][k]); + ofs << " " << std::setw(10) << val; + } + ofs << std::endl; + } + } + else + { + // Real 2m x 2m matrix in spin basis; this is the + // layout parsed by read_occup_m() on restart. + for (int m0 = 0; m0 < m; ++m0) + { + for (int ipol0 = 0; ipol0 < npol; ++ipol0) + { + for (int m1 = 0; m1 < m; ++m1) { - for (int ipol1 = 0; ipol1 < npol; ipol1++) + for (int ipol1 = 0; ipol1 < npol; ++ipol1) { - int m1_all = m1 + (2 * l + 1) * ipol1; - ofs << std::setw(12) << std::setprecision(8) << std::fixed - << dftu.occmat().get(iat, l, 0, m0_all, m1_all); + const int k = m0 * m + m1; + double val = 0.0; + if (ipol0 == 0 && ipol1 == 0) + { + val = 0.5 * (blocks[0][k] + blocks[3][k]); + } + else if (ipol0 == 1 && ipol1 == 1) + { + val = 0.5 * (blocks[0][k] - blocks[3][k]); + } + else + { + val = 0.5 * blocks[1][k]; + } + ofs << " " << std::setw(10) + << std::setprecision(8) + << std::fixed << val; } } ofs << std::endl; diff --git a/source/source_pw/module_pwdft/dftu_base_io.h b/source/source_pw/module_pwdft/dftu_base_io.h index efa721d74b6..8b24b39eede 100644 --- a/source/source_pw/module_pwdft/dftu_base_io.h +++ b/source/source_pw/module_pwdft/dftu_base_io.h @@ -3,6 +3,7 @@ #include "source_base/matrix.h" #include "source_estate/occ_matrix.h" +#include "source_pw/module_pwdft/dftu_types.h" #include #include @@ -17,19 +18,97 @@ namespace DFTU_BASE /// nested occupation-matrix type used by DFT+U: occ_mat[iat][l][spin](m0, m1) using OccMatData = std::vector>>; +/// Tell whether the given ionic step produces a numbered file. +/// +/// @param istep ionic-step index, starting from 0 +/// @param cfg frequency configuration +/// @return true when out_freq_ion > 0 and istep is an output ionic step +bool is_ion_step_output_step(int istep, const OccmatOutputCfg& cfg); + +/// Tell whether the current electronic step is recorded inside the numbered file. +/// +/// @param iter electronic-iteration index, starting from 1 +/// @param conv_esolver whether the electronic SCF is converged at this step +/// @param cfg frequency configuration +/// @return true when iter % out_freq_elec == 0, iter == scf_nmax, or converged +bool is_elec_snapshot_trigger(int iter, + bool conv_esolver, + const OccmatOutputCfg& cfg); + +/// Build the name of the per-ionic-step occupation-matrix file. +/// +/// @param out_dir output directory (including the trailing separator) +/// @param istep ionic-step index, starting from 0 +/// @return occ_matg{istep+1}.txt +std::string gen_ion_step_occ_mat_filename(const std::string& out_dir, int istep); + +/// Return the full path of the first candidate file that exists in @p dir. +/// Returns an empty string if none of the candidates exist. +/// Only the calling process probes the filesystem; callers must arrange +/// MPI broadcast of the result if other ranks need it. +std::string find_first_existing_file(const std::string& dir, + const std::vector& candidates); + +/// Append one electronic-step section to the per-ionic-step file. +/// +/// The section records the electronic-step index, the configured charge- +/// density convergence threshold (scf_thr) and the actual residual (drho) +/// of the current electronic step, followed by the occupation matrices and +/// per-atom magnetism. When occmat_ready is false, an "N/A" placeholder is +/// recorded instead of the matrix body (the PW path has no matrix at +/// istep 0 / iter 1 unless it was loaded from file). The caller is +/// responsible for truncating the file at the first electronic step of the +/// ionic step (see prepare_ion_step_file()) and for invoking this function +/// only after drho of the current electronic step is computed. +/// +/// @param dftu DFT+U object holding the occupation matrices +/// @param ucell unit cell +/// @param global_out_dir output directory (including the trailing separator) +/// @param nspin number of spin components (1, 2 or 4) +/// @param npol number of polarizations +/// @param istep ionic-step index, starting from 0 +/// @param iter electronic-iteration index, starting from 1 +/// @param conv_esolver whether the electronic SCF is converged at this step +/// @param occmat_ready whether the occupation matrix of this step exists; +/// false records an "N/A" placeholder instead of the matrix body +/// @param scf_thr configured charge-density convergence threshold +/// @param drho actual charge-density residual of the current electronic step +/// @param cfg frequency configuration +/// @param soc_layout storage layout of the nspin == 4 occupation matrix +void append_ion_step_snapshot(const Plus_U_Base& dftu, + const UnitCell& ucell, + const std::string& global_out_dir, + int nspin, + int npol, + int istep, + int iter, + bool conv_esolver, + bool occmat_ready, + double scf_thr, + double drho, + const OccmatOutputCfg& cfg, + OccmatSocLayout soc_layout); + /// 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 -/// occ_mat_ctrl and init_chg. +/// init_occ_mat and init_chg. +/// +/// @param soc_layout storage layout of the nspin == 4 occupation matrix; +/// the PW path stores 4 contiguous Pauli blocks (b0, b1, b2, b3) and +/// must reconstruct them from the (n_uu, Re(n_ud), Im(n_ud), n_dd) +/// blocks written to the file; the LCAO path stores a real 2m x 2m +/// spin-basis matrix and discards Im(n_ud). void read_occup_m(const UnitCell& ucell, OccupationMatrix& occ, const std::vector& l_channel, - const int occ_mat_ctrl, + const int init_occ_mat, const std::string& fn, const std::string& init_chg, int nspin, - int npol); + int npol, + OccmatSocLayout soc_layout); /// Broadcast the local occupation number matrices from rank 0 to all ranks. /// @@ -40,30 +119,92 @@ void local_occup_bcast(const UnitCell& ucell, int nspin, int npol); +/// Create (or truncate) the per-ionic-step file occ_matg{istep+1}.txt and +/// write its header (rank 0 only). +/// +/// Must be called once at the first electronic step (iter == 1) of an output +/// ionic step, before any append_ion_step_snapshot() section. Truncating here +/// also guarantees that a rerun in the same output directory can not append +/// snapshots of a previous calculation. +void prepare_ion_step_file(const std::string& global_out_dir, + const int istep, + const OccmatOutputCfg& cfg); + /// Output DFT+U information (Hubbard U/J, local occupation matrices) to the -/// running log and, when out_chg is set, to the dm_onsite.txt file. +/// running log. +/// +/// When cfg.out_occ_mat is true, cfg.out_freq_ion is positive and istep is an +/// output ionic step, the per-ionic-step file occ_matg{istep+1}.txt is +/// created: at the first electronic step (iter == 1) it is truncated and +/// initialized with a provenance header, then append_ion_step_snapshot() +/// appends one section per recorded electronic step. +/// +/// Note: occ_mat.txt is NOT written here. It records the actual +/// charge-density residual drho, which is only known after the electronic +/// solve, so write_latest_occmat() writes it from the iter_finish stage. /// /// Extracted from Plus_U_Base::output as a free function so that IO logic is /// decoupled from the Plus_U_Base class. The function only reads the /// Plus_U_Base state via public accessors; no friend declaration needed. void output(const Plus_U_Base& dftu, const UnitCell& ucell, - bool out_chg, const std::string& global_out_dir, int nspin, - int npol); + int npol, + int istep, + int iter, + const OccmatOutputCfg& cfg, + OccmatSocLayout soc_layout); + +/// Overwrite occ_mat.txt with the occupation matrix of the current +/// electronic step (rank 0 only). +/// +/// occ_mat.txt is a single-section snapshot file: the same provenance +/// header and compact layout as occ_matg{#}.txt, whose section header +/// additionally carries the configured scf_thr and the actual drho of +/// this step. Must be called at the iter_finish stage, after drho is +/// computed; it is the entry file of init_chg=file and NSCF restarts. +/// No-op when cfg.out_occ_mat is false. +/// +/// @param dftu DFT+U object holding the occupation matrices +/// @param ucell unit cell +/// @param global_out_dir output directory (including the trailing separator) +/// @param nspin number of spin components (1, 2 or 4) +/// @param npol number of polarizations +/// @param istep ionic-step index, starting from 0 +/// @param iter electronic-iteration index, starting from 1 +/// @param scf_thr configured charge-density convergence threshold +/// @param drho actual charge-density residual of the current electronic step +/// @param cfg output configuration; only cfg.out_occ_mat is consulted here +/// @param soc_layout storage layout of the nspin == 4 occupation matrix +void write_latest_occmat(const Plus_U_Base& dftu, + const UnitCell& ucell, + const std::string& global_out_dir, + int nspin, + int npol, + int istep, + int iter, + double scf_thr, + double drho, + const OccmatOutputCfg& cfg, + OccmatSocLayout soc_layout); /// Write local occupation matrices to the given stream. /// -/// Extracted from Plus_U_Base::write_occup_m. When diag is true, eigenvalues -/// and magnetism are also printed; otherwise only raw matrix elements. +/// When diag is true, eigenvalues and magnetism are also printed; otherwise +/// only raw matrix elements. fmt selects between the legacy token layout +/// (used for running log) and the readable snapshot layout of +/// occ_mat.txt and occ_matg{#}.txt. soc_layout tells how the nspin == 4 +/// storage is arranged (PW Pauli blocks or LCAO real spin-basis matrix). /// Caller is responsible for opening/closing the stream. void write_occup_m(const Plus_U_Base& dftu, const UnitCell& ucell, - std::ofstream& ofs, + std::ostream& ofs, bool diag, int nspin, - int npol); + int npol, + OccmatTextFormat fmt, + OccmatSocLayout soc_layout); } // namespace DFTU_BASE diff --git a/source/source_pw/module_pwdft/dftu_types.h b/source/source_pw/module_pwdft/dftu_types.h new file mode 100644 index 00000000000..f409e5b884b --- /dev/null +++ b/source/source_pw/module_pwdft/dftu_types.h @@ -0,0 +1,46 @@ +#ifndef DFTU_TYPES_H +#define DFTU_TYPES_H + +namespace DFTU_BASE +{ + +/// Frequency configuration for numbered occupation-matrix output. +/// +/// Each output ionic step (out_freq_ion) owns one file occ_matg{#}.txt. +/// Within that file, one section is appended at every electronic step +/// selected by out_freq_elec, scf_nmax or convergence. +struct OccmatOutputCfg +{ + int out_freq_ion; ///< ionic-step interval; 0 disables numbered files + int out_freq_elec; ///< electronic-iteration interval recorded inside a numbered file + int scf_nmax; ///< maximum number of electronic iterations + bool out_occ_mat; ///< master switch of occupation-matrix output (out_occ_mat INPUT parameter) + /// DFT+U calculation mode (INPUT parameter dft_plus_u). The occupation- + /// matrix IO is meaningful only when dft_plus_u > 0; otherwise the + /// occupation matrix is never computed and the writers must early-out + /// to avoid touching an empty l_channel vector. This honours the + /// documented contract that out_occ_mat only takes effect for DFT+U. + int dft_plus_u; +}; + +/// Text format used by write_occup_m(). +enum OccmatTextFormat +{ + OCMAT_FMT_LEGACY, ///< whitespace-token layout parsed by read_occup_m() + OCMAT_FMT_READABLE ///< compact human-readable layout for snapshot files +}; + +/// Storage layout of the real (2l+1)*npol matrix when nspin == 4. +enum OccmatSocLayout +{ + /// 4 contiguous Pauli-component blocks [charge, sigma_x, sigma_y, sigma_z], + /// each (2l+1)x(2l+1); used by the PW path + SOC_LAYOUT_PAULI, + /// real symmetric matrix in spin basis (m + spin-polarization index); + /// imaginary parts are not stored; used by the LCAO path + SOC_LAYOUT_SPIN_BASIS_REAL +}; + +} // namespace DFTU_BASE + +#endif // DFTU_TYPES_H diff --git a/source/source_pw/module_pwdft/exx_helper.cpp b/source/source_pw/module_pwdft/exx_helper.cpp index 3f117482540..b183fc0e077 100644 --- a/source/source_pw/module_pwdft/exx_helper.cpp +++ b/source/source_pw/module_pwdft/exx_helper.cpp @@ -58,7 +58,7 @@ void Exx_Helper::before_scf(void* p_hamilt, void* psi, const Input_pa template bool Exx_Helper::iter_finish(void* p_elec, Charge* p_charge, void* psi, UnitCell& ucell, const Input_para& inp, - bool& conv_esolver, int& iter) + bool& conv_esolver, int iter) { /// Return if EXX is not enabled if (op_exx == nullptr) @@ -125,7 +125,7 @@ double Exx_Helper::cal_exx_energy(void* psi_) } template -bool Exx_Helper::exx_after_converge(int &iter, bool ene_conv) +bool Exx_Helper::exx_after_converge(int iter, bool ene_conv) { if (op_exx->first_iter) { @@ -150,7 +150,13 @@ bool Exx_Helper::exx_after_converge(int &iter, bool ene_conv) return true; } GlobalV::ofs_running << "Updating EXX and rerun SCF" << std::endl; - iter = 0; + // SCF rerun is now signaled by returning false (not converged). The + // caller (ESolver_KS_PW::iter_finish) detects "SCF was converged but + // EXX says rerun" and sets ESolver_KS::scf_rerun_, which makes + // ESolver_KS::runner restart the loop at iter=1. The old trick of + // writing iter=0 here is removed: it crashed DFT+U occupation-matrix + // writers (append_ion_step_snapshot validates iter>=1) and broke the + // 1-based iteration invariant for every downstream reader of iter. return false; } diff --git a/source/source_pw/module_pwdft/exx_helper.h b/source/source_pw/module_pwdft/exx_helper.h index 2adb945b8b6..31b25eeeadf 100644 --- a/source/source_pw/module_pwdft/exx_helper.h +++ b/source/source_pw/module_pwdft/exx_helper.h @@ -26,7 +26,7 @@ struct Exx_Helper : public Exx_HelperBase bool iter_finish(void* p_elec, Charge* p_charge, void* psi, UnitCell& ucell, const Input_para& inp, - bool& conv_esolver, int& iter) override; + bool& conv_esolver, int iter) override; void set_firstiter(bool flag = true) override { first_iter = flag; } void set_wg(const ModuleBase::matrix *wg_) override { wg = wg_; } @@ -40,7 +40,7 @@ struct Exx_Helper : public Exx_HelperBase op_exx->set_wg(wg); } - bool exx_after_converge(int &iter, bool ene_conv) override; + bool exx_after_converge(int iter, bool ene_conv) override; double cal_exx_energy(void* psi_) override; diff --git a/source/source_pw/module_pwdft/exx_helper_base.h b/source/source_pw/module_pwdft/exx_helper_base.h index 7ebab443948..12bd8d7aeea 100644 --- a/source/source_pw/module_pwdft/exx_helper_base.h +++ b/source/source_pw/module_pwdft/exx_helper_base.h @@ -20,7 +20,7 @@ class Exx_HelperBase virtual bool iter_finish(void* p_elec, Charge* p_charge, void* psi, UnitCell& ucell, const Input_para& inp, - bool& conv_esolver, int& iter) = 0; + bool& conv_esolver, int iter) = 0; virtual void set_firstiter(bool flag = true) = 0; virtual void set_wg(const ModuleBase::matrix* wg) = 0; @@ -29,7 +29,13 @@ class Exx_HelperBase virtual void set_op() = 0; - virtual bool exx_after_converge(int& iter, bool ene_conv) = 0; + /// @brief Decide whether EXX has converged and SCF should rerun. + /// @note iter is taken by value: EXX no longer mutates the caller's + /// iteration counter. SCF rerun is requested via the return + /// value (false = not converged, rerun SCF) and the caller + /// (ESolver_KS_PW) sets ESolver_KS::scf_rerun_ to trigger the + /// restart at iter=1 in ESolver_KS::runner. + virtual bool exx_after_converge(int iter, bool ene_conv) = 0; virtual double cal_exx_energy(void* psi) = 0; diff --git a/source/source_pw/module_pwdft/setup_dftu_pw.cpp b/source/source_pw/module_pwdft/setup_dftu_pw.cpp index 98cafa75f28..6e0332eac87 100644 --- a/source/source_pw/module_pwdft/setup_dftu_pw.cpp +++ b/source/source_pw/module_pwdft/setup_dftu_pw.cpp @@ -2,7 +2,7 @@ #include "source_pw/module_pwdft/dftu_base.h" // mohan add 2025-11-06 #include "source_pw/module_pwdft/dftu_base_io.h" // mohan add 2025-11-08 #include "source_pw/module_pwdft/dftu_pw.h" -#include "source_io/module_parameter/parameter.h" +#include "source_pw/module_pwdft/dftu_pw_tools.h" namespace DFTU_BASE { @@ -14,29 +14,77 @@ void iter_init_dftu_pw(const int iter, const ModuleBase::matrix& wg, const UnitCell& ucell, Charge_Mixing* p_chgmix, - const int* isk) + const std::string& global_out_dir, + const OccmatOutputCfg& occmat_cfg, + const int* isk, + const int nspin, + const int kpar, + const int npol) { - if (!p_chgmix || !PARAM.inp.dft_plus_u) + if (!p_chgmix || occmat_cfg.dft_plus_u <= 0) { return; } - if (iter == 1 && istep == 0) + // Prepare the per-ionic-step file at the first electronic step, before + // the first occupation matrix exists. output() is skipped below at + // istep 0 / iter 1; without this call the g1 file would miss its header + // and could inherit stale sections from a previous run. + if (iter == 1 && occmat_cfg.out_occ_mat + && DFTU_BASE::is_ion_step_output_step(istep, occmat_cfg)) { - return; + DFTU_BASE::prepare_ion_step_file(global_out_dir, + istep, occmat_cfg); } - if (dftu.get_occ_mat_ctrl() != 2) + const bool first_scf_step = (iter == 1 && istep == 0); + + // Build the +U potential from the file-loaded occupation matrix when: + // - init_occ_mat == 2: the matrix is fixed for the whole run, so + // rebuild pot_onsite and E_plusU from it at every electronic step + // (the matrix is never reprojected from psi or mixed); + // - init_occ_mat == 1 at the first SCF step: the file matrix seeds + // the Hamiltonian so the first diagonalization includes the +U + // correction. On later steps cal_occ_pw reprojects the matrix + // from psi (which now carries the +U influence) and rebuilds the + // potential, matching the documented "first step from file, then + // update" semantics. Without this the file matrix would never + // enter the Hamiltonian (default mixing_dftu == 0). + if (dftu.get_init_occ_mat() == 2 + || (first_scf_step && dftu.get_init_occ_mat() == 1)) { - DFTU_BASE::cal_occ_pw(psi, wg, ucell, p_chgmix, isk, PARAM.inp.kpar, - PARAM.inp.nspin, dftu.get_device(), + pw::compute_pot_uterm_and_energy(ucell, + nspin, + dftu.get_u_current_vec(), + dftu.get_l_channel_vec(), + dftu.get_uterm_mat_index(), + dftu.occmat(), + dftu.get_uterm_mat(), + dftu.energy_ref()); + } + else if (!first_scf_step) + { + // No occupation matrix exists before the first projection, so the + // matrix can only be accumulated from psi on later electronic steps. + DFTU_BASE::cal_occ_pw(psi, wg, ucell, p_chgmix, isk, kpar, + nspin, dftu.get_device(), dftu.get_l_channel_vec(), dftu.get_u_current_vec(), dftu.get_uterm_mat_index(), dftu.occmat(), dftu.has_occ_mixer() ? &dftu.occ_mixer() : nullptr, dftu.get_uterm_mat(), dftu.energy_ref()); } - DFTU_BASE::output(dftu, ucell, PARAM.inp.out_chg[0], PARAM.globalv.global_out_dir, PARAM.inp.nspin, PARAM.globalv.npol); + + // Keep the historical output cadence: output() is skipped at the very + // first electronic step. + if (first_scf_step) + { + return; + } + + DFTU_BASE::output(dftu, ucell, global_out_dir, + nspin, npol, istep, iter, occmat_cfg, + DFTU_BASE::SOC_LAYOUT_PAULI); } } diff --git a/source/source_pw/module_pwdft/setup_dftu_pw.h b/source/source_pw/module_pwdft/setup_dftu_pw.h index fa6c43bc74a..05c448e2164 100644 --- a/source/source_pw/module_pwdft/setup_dftu_pw.h +++ b/source/source_pw/module_pwdft/setup_dftu_pw.h @@ -5,12 +5,16 @@ #include "source_base/matrix.h" #include "source_estate/module_charge/chg_mix.h" +#include + struct Input_para; class Plus_U_Base; // mohan add 2025-11-06 namespace DFTU_BASE { +struct OccmatOutputCfg; + void iter_init_dftu_pw(const int iter, const int istep, Plus_U_Base& dftu, // mohan add 2025-11-06 @@ -18,7 +22,12 @@ void iter_init_dftu_pw(const int iter, const ModuleBase::matrix& wg, const UnitCell& ucell, Charge_Mixing* p_chgmix, - const int* isk); + const std::string& global_out_dir, + const OccmatOutputCfg& occmat_cfg, + const int* isk, + const int nspin, + const int kpar, + const int npol); } diff --git a/source/source_pw/module_pwdft/setup_pot.cpp b/source/source_pw/module_pwdft/setup_pot.cpp index 44da50a44e0..eec2cf13ad6 100644 --- a/source/source_pw/module_pwdft/setup_pot.cpp +++ b/source/source_pw/module_pwdft/setup_pot.cpp @@ -130,8 +130,9 @@ void pw::setup_pot(const int istep, inp.device, PARAM.globalv.hubbard_u, PARAM.globalv.uramping, - inp.occ_mat_ctrl, - inp.mixing_dftu); + inp.init_occ_mat, + inp.mixing_dftu, + DFTU_BASE::SOC_LAYOUT_PAULI); } return; diff --git a/source/source_pw/module_pwdft/test/CMakeLists.txt b/source/source_pw/module_pwdft/test/CMakeLists.txt index 93b877b7682..ad7674a93e8 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_pw_tools.cpp ../yukawa_screening.cpp ../../../source_estate/occ_matrix.cpp ../../../source_estate/occ_mixer.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 cb2c3dc3a0d..aca4b5f0e45 100644 --- a/source/source_pw/module_pwdft/test/dftu_base_test.cpp +++ b/source/source_pw/module_pwdft/test/dftu_base_test.cpp @@ -9,12 +9,22 @@ ***********************************************/ #include "source_pw/module_pwdft/dftu_base.h" +#include "source_pw/module_pwdft/dftu_base_io.h" +#include "source_pw/module_pwdft/dftu_pw_tools.h" #include "source_cell/atom_spec.h" #include "source_cell/unitcell.h" #include "gtest/gtest.h" +#ifdef __MPI +#include +#endif + +#include +#include +#include +#include #include #include @@ -91,8 +101,9 @@ class DFTUBaseTest : public testing::Test "cpu", // device hubbard_u, 0.0, // uramping - 0, // occ_mat_ctrl - 0); // mixing_dftu + 0, // init_occ_mat + 0, // mixing_dftu + DFTU_BASE::SOC_LAYOUT_PAULI); } }; @@ -378,3 +389,659 @@ TEST_F(OccMatRoundtripTest, Nspin4_PauliBlocks) for (int j = 0; j < sizes[i]; j++) EXPECT_DOUBLE_EQ(occ_mat[i].data[j], static_cast(i * 1000 + j + 1)); } + +/// append_ion_step_snapshot must record an "N/A" placeholder instead of a +/// silent zero matrix when the occupation matrix does not exist yet (the PW +/// istep 0 / iter 1 case), and the real matrix body when it does. +TEST_F(DFTUBaseTest, AppendSnapshotNAPlaceholderAndReady) +{ + Plus_U_Base dftu; + this->init_dftu(dftu, false); + + const DFTU_BASE::OccmatOutputCfg cfg = {1, 1, 5, true, 1}; + const std::string out_dir = testing::TempDir(); + + // fresh run: no occupation matrix has been computed or loaded + ASSERT_FALSE(dftu.is_occmat_ready()); + DFTU_BASE::append_ion_step_snapshot(dftu, + ucell, + out_dir, + 2, // nspin + 1, // npol + 0, // istep -> occ_matg1.txt + 1, // iter + false, + false, // occmat_ready + 1e-6, + 0.5, + cfg, + DFTU_BASE::SOC_LAYOUT_PAULI); + + const std::string fn1 = out_dir + "occ_matg1.txt"; + std::ifstream ifs(fn1.c_str()); + ASSERT_TRUE(ifs.is_open()); + std::stringstream ss; + ss << ifs.rdbuf(); + ifs.close(); + const std::string content = ss.str(); + EXPECT_NE(content.find("# Electronic step 1"), std::string::npos); + EXPECT_NE(content.find("# scf_thr 1.00000000e-06"), std::string::npos); + EXPECT_NE(content.find("# drho 5.00000000e-01"), std::string::npos); + EXPECT_NE(content.find("\n N/A\n"), std::string::npos); + EXPECT_EQ(content.find("Atom"), std::string::npos); + + // ready case: the matrix body is written for the atom (fresh file g2) + DFTU_BASE::append_ion_step_snapshot(dftu, + ucell, + out_dir, + 2, // nspin + 1, // npol + 1, // istep -> occ_matg2.txt + 2, // iter + false, + true, // occmat_ready + 1e-6, + 0.5, + cfg, + DFTU_BASE::SOC_LAYOUT_PAULI); + + const std::string fn2 = out_dir + "occ_matg2.txt"; + std::ifstream ifs2(fn2.c_str()); + ASSERT_TRUE(ifs2.is_open()); + std::stringstream ss2; + ss2 << ifs2.rdbuf(); + ifs2.close(); + const std::string content2 = ss2.str(); + EXPECT_NE(content2.find("Fe Atom 1 L 2"), std::string::npos); + EXPECT_EQ(content2.find("N/A"), std::string::npos); + + // Clean up so subsequent tests checking "file not written" are not + // fooled by leftovers. std::remove is Rule 17 compliant. + std::remove((out_dir + "occ_matg1.txt").c_str()); + std::remove((out_dir + "occ_matg2.txt").c_str()); +} + +/// out_occ_mat = false must suppress both the numbered snapshot file and the +/// latest occ_mat.txt, even when the frequency gates would trigger. +TEST_F(DFTUBaseTest, OccMatSwitchDisabledWritesNothing) +{ + Plus_U_Base dftu; + this->init_dftu(dftu, false); + + const DFTU_BASE::OccmatOutputCfg cfg = {1, 1, 5, false, 1}; + const std::string out_dir = testing::TempDir(); + + // Remove stale files from previous test binary runs so the + // "file not written" assertions are meaningful. + std::remove((out_dir + "occ_matg1.txt").c_str()); + std::remove((out_dir + "occ_mat.txt").c_str()); + + DFTU_BASE::append_ion_step_snapshot(dftu, + ucell, + out_dir, + 2, // nspin + 1, // npol + 0, // istep (an output ionic step) + 1, // iter + false, + true, // occmat_ready + 1e-6, + 0.5, + cfg, + DFTU_BASE::SOC_LAYOUT_PAULI); + + const std::string fn_g1 = out_dir + "occ_matg1.txt"; + std::ifstream ifs(fn_g1.c_str()); + EXPECT_FALSE(ifs.is_open()); + + DFTU_BASE::write_latest_occmat(dftu, + ucell, + out_dir, + 2, // nspin + 1, // npol + 0, // istep + 2, // iter + 1e-6, + 0.5, + cfg, + DFTU_BASE::SOC_LAYOUT_PAULI); + + const std::string fn_latest = out_dir + "occ_mat.txt"; + std::ifstream ifs_latest(fn_latest.c_str()); + EXPECT_FALSE(ifs_latest.is_open()); +} + +/// dft_plus_u = 0 (no DFT+U) must suppress both the numbered snapshot file +/// and the latest occ_mat.txt, even when out_occ_mat = true (the default). +/// This is the regression test for the abacuslite segfault: a default- +/// constructed Plus_U_Base has an empty l_channel vector, so without this +/// guard write_occup_m() dereferences a null l_channel.data() inside +/// has_l_channel(). The guard also honours the documented contract that +/// out_occ_mat only takes effect for DFT+U calculations (dft_plus_u > 0). +TEST_F(DFTUBaseTest, DftPlusUDisabledWritesNothing) +{ + Plus_U_Base dftu; // default-constructed: l_channel is empty, mirroring + // the LCAO/PW esolver path when dft_plus_u == 0 + // Deliberately skip init_dftu(): the bug is that the IO functions must + // not even reach has_l_channel() when dft_plus_u == 0. + + const DFTU_BASE::OccmatOutputCfg cfg = {1, 1, 5, true, 0}; + const std::string out_dir = testing::TempDir() + "guard_test/"; + + DFTU_BASE::append_ion_step_snapshot(dftu, + ucell, + out_dir, + 2, // nspin + 1, // npol + 0, // istep (an output ionic step) + 1, // iter + false, + true, // occmat_ready + 1e-6, + 0.5, + cfg, + DFTU_BASE::SOC_LAYOUT_PAULI); + + const std::string fn_g1 = out_dir + "occ_matg1.txt"; + std::ifstream ifs(fn_g1.c_str()); + EXPECT_FALSE(ifs.is_open()); + + DFTU_BASE::write_latest_occmat(dftu, + ucell, + out_dir, + 2, // nspin + 1, // npol + 0, // istep + 2, // iter + 1e-6, + 0.5, + cfg, + DFTU_BASE::SOC_LAYOUT_PAULI); + + const std::string fn_latest = out_dir + "occ_mat.txt"; + std::ifstream ifs_latest(fn_latest.c_str()); + EXPECT_FALSE(ifs_latest.is_open()); +} + +/// find_first_existing_file must return the first candidate that exists, +/// in declaration order, and an empty string when none exist. +TEST(FindFirstExistingFileTest, ReturnsFirstExistingCandidate) +{ + // Use the gtest-managed temp dir for writable fixtures, and a + // non-existent subdir for the "no candidate" case. The test never + // creates or deletes directories itself (AGENTS.md rule 17). + const std::string td = testing::TempDir(); + + // No candidates exist -> empty string. All three names are unique to + // this test, so they should not be present in td. + { + const std::vector candidates = {"ffe_absent_1.txt", + "ffe_absent_2.txt", + "ffe_absent_3.txt"}; + EXPECT_TRUE(DFTU_BASE::find_first_existing_file(td, candidates).empty()); + } + + // Only the second candidate exists -> it is returned even though the + // first comes earlier in the list. + { + const std::string fn = td + "ffe_second_only.txt"; + std::ofstream(fn).close(); + const std::vector candidates = {"ffe_absent_1.txt", + "ffe_second_only.txt", + "ffe_absent_3.txt"}; + EXPECT_EQ(DFTU_BASE::find_first_existing_file(td, candidates), fn); + } + + // The first candidate exists -> it takes precedence over the second. + { + const std::string fn = td + "ffe_first_wins.txt"; + std::ofstream(fn).close(); + const std::vector candidates = {"ffe_first_wins.txt", + "ffe_second_only.txt", + "ffe_absent_3.txt"}; + EXPECT_EQ(DFTU_BASE::find_first_existing_file(td, candidates), fn); + } +} + +/// For init_occ_mat=2 the occupation-matrix file is read exactly once. +/// Later init_base() calls (one per ionic step in a relax run) must keep +/// the in-memory matrix instead of reading the file again, so pointing +/// the second call at a non-existent readin dir must not abort the run. +TEST_F(DFTUBaseTest, InitBaseReadsOccMatFileOnlyOnce) +{ + // Use the gtest-managed temporary directory so the test never creates + // or deletes directories itself (see AGENTS.md rule 17). + const std::string dir = testing::TempDir(); + const std::string fn = dir + "occ_mat.txt"; + + // One Fe atom, L=2, two spin channels of 5x5, filled with distinct + // constants so a re-read would be easy to distinguish from a + // preserved in-memory matrix. + { + std::ofstream ofs(fn); + ASSERT_TRUE(ofs.is_open()); + ofs << "# compact test fixture\n"; + ofs << " Fe Atom 1 L 2 mag 0.0\n"; + ofs << " spin 1 nelec 0.5\n"; + for (int m0 = 0; m0 < 5; ++m0) + { + for (int m1 = 0; m1 < 5; ++m1) + { + ofs << " 0.1"; + } + ofs << "\n"; + } + ofs << " spin 2 nelec 0.5\n"; + for (int m0 = 0; m0 < 5; ++m0) + { + for (int m1 = 0; m1 < 5; ++m1) + { + ofs << " 0.2"; + } + ofs << "\n"; + } + } + + Plus_U_Base dftu; + const std::vector l_channel = {2}; + const std::vector hubbard_u = {0.0}; + auto call_init = [&](const std::string& readin_dir) + { + dftu.init_base(ucell, + 1, // npol + 2, // nspin + l_channel, + false, // yukawa_potential + 0.5, // yukawa_lambda + readin_dir, // global_readin_dir + "", // global_out_dir + "none", // init_chg + "cpu", // device + hubbard_u, + 0.0, // uramping + 2, // init_occ_mat + 0, // mixing_dftu + DFTU_BASE::SOC_LAYOUT_PAULI); + }; + + // First ionic step: the file is read. + call_init(dir); + ASSERT_TRUE(dftu.is_occmat_ready()); + EXPECT_NEAR(dftu.occmat().get(0, 2, 0, 0, 0), 0.1, 1e-12); + EXPECT_NEAR(dftu.occmat().get(0, 2, 1, 4, 4), 0.2, 1e-12); + + // Simulate a later ionic step: the in-memory matrix is preserved and + // the file is not consulted again, so a non-existent readin dir is fine. + call_init(dir + "does_not_exist/"); + EXPECT_TRUE(dftu.is_occmat_ready()); + EXPECT_NEAR(dftu.occmat().get(0, 2, 0, 0, 0), 0.1, 1e-12); + EXPECT_NEAR(dftu.occmat().get(0, 2, 1, 4, 4), 0.2, 1e-12); +} + +/// Reading an nspin=4 occupation-matrix file with SOC_LAYOUT_PAULI must +/// reconstruct the 4 contiguous Pauli blocks [b0, b1, b2, b3] in the 2m x 2m +/// flat buffer. The bug fixed in commit 1787365b3 was that the Im(n_ud) +/// block ("spin 12 im") was discarded, so b2 came out zero/garbage and the +/// subsequent SCF diverged. This test writes a fixture with a non-zero Im +/// block and verifies all 4 blocks at the correct flat-buffer offsets. +TEST_F(DFTUBaseTest, InitBaseReadsOccMatSocPauliRoundtrip) +{ + // Use the gtest-managed temporary directory (AGENTS.md rule 17). + const std::string dir = testing::TempDir(); + const std::string fn = dir + "occ_mat.txt"; + + // One Fe atom, L=2 -> nm=5. The 2m x 2m flat buffer (100 elements) + // holds 4 contiguous m^2 blocks: [b0, b1, b2, b3]. The file stores + // (n_uu, Re(n_ud), Im(n_ud), n_dd); the reader must reconstruct + // b0 = n_uu + n_dd + // b1 = 2 * Re(n_ud) + // b2 = 2 * Im(n_ud) <- discarded by the buggy reader + // b3 = n_uu - n_dd + const int nm = 5; + const int m2 = nm * nm; + std::vector uu(m2), re(m2), im(m2), dd(m2); + for (int k = 0; k < m2; ++k) + { + uu[k] = 1.0 + 0.01 * k; + re[k] = 0.1 + 0.01 * k; + im[k] = 0.01 + 0.001 * k; // non-zero -- the bug discarded this block + dd[k] = 0.5 + 0.01 * k; + } + + { + 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"; + for (int m0 = 0; m0 < nm; ++m0) + { + for (int m1 = 0; m1 < nm; ++m1) + ofs << " " << uu[m0 * nm + m1]; + ofs << "\n"; + } + ofs << " spin 12 re\n"; + for (int m0 = 0; m0 < nm; ++m0) + { + for (int m1 = 0; m1 < nm; ++m1) + ofs << " " << re[m0 * nm + m1]; + ofs << "\n"; + } + ofs << " spin 12 im\n"; + for (int m0 = 0; m0 < nm; ++m0) + { + for (int m1 = 0; m1 < nm; ++m1) + ofs << " " << im[m0 * nm + m1]; + ofs << "\n"; + } + ofs << " spin 2 nelec 0.5\n"; + for (int m0 = 0; m0 < nm; ++m0) + { + for (int m1 = 0; m1 < nm; ++m1) + ofs << " " << dd[m0 * nm + m1]; + ofs << "\n"; + } + } + + Plus_U_Base dftu; + const std::vector l_channel = {2}; + const std::vector hubbard_u = {0.0}; + dftu.init_base(ucell, + 2, // npol (nspin=4 requires npol=2) + 4, // nspin + l_channel, + false, // yukawa_potential + 0.5, // yukawa_lambda + dir, // global_readin_dir + "", // global_out_dir + "none", // init_chg + "cpu", // device + hubbard_u, + 0.0, // uramping + 2, // init_occ_mat + 0, // mixing_dftu + DFTU_BASE::SOC_LAYOUT_PAULI); + + ASSERT_TRUE(dftu.is_occmat_ready()); + + // The 2m x 2m flat buffer holds 4 contiguous m^2 blocks: [b0, b1, b2, b3]. + const ModuleBase::matrix& occ0 = dftu.occmat().mat(0, 2, 0); + ASSERT_EQ(occ0.nr * occ0.nc, 4 * m2); + + for (int k = 0; k < m2; ++k) + { + EXPECT_NEAR(occ0.c[0 * m2 + k], uu[k] + dd[k], 1e-12) // b0 + << "Pauli block b0 mismatch at k=" << k; + EXPECT_NEAR(occ0.c[1 * m2 + k], 2.0 * re[k], 1e-12) // b1 + << "Pauli block b1 mismatch at k=" << k; + EXPECT_NEAR(occ0.c[2 * m2 + k], 2.0 * im[k], 1e-12) // b2 -- the bug + << "Pauli block b2 mismatch at k=" << k; + EXPECT_NEAR(occ0.c[3 * m2 + k], uu[k] - dd[k], 1e-12) // b3 + << "Pauli block b3 mismatch at k=" << k; + } +} + +/// 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 +/// E_plusU stayed zero and the file matrix was overwritten by the +/// psi-reprojected matrix at iter=2 without ever entering the Hamiltonian. +/// This test verifies the mechanism: after init_base reads the file, +/// compute_pot_uterm_and_energy must produce a non-zero potential and +/// energy from the file-loaded matrix. +TEST_F(DFTUBaseTest, InitOccMat1BuildsPotentialFromFile) +{ + const std::string dir = testing::TempDir(); + const std::string fn = dir + "occ_mat.txt"; + + // One Fe atom, L=2, nspin=2. All 25 elements per spin set to a + // constant so expected pot_onsite and energy are easy to compute. + const double occ_up = 0.1; + const double occ_dn = 0.2; + { + std::ofstream ofs(fn); + ASSERT_TRUE(ofs.is_open()); + ofs << "# compact test fixture\n"; + ofs << " Fe Atom 1 L 2 mag 0.0\n"; + ofs << " spin 1 nelec 0.5\n"; + for (int m0 = 0; m0 < 5; ++m0) + { + for (int m1 = 0; m1 < 5; ++m1) { ofs << " " << occ_up; } + ofs << "\n"; + } + ofs << " spin 2 nelec 0.5\n"; + for (int m0 = 0; m0 < 5; ++m0) + { + for (int m1 = 0; m1 < 5; ++m1) { ofs << " " << occ_dn; } + ofs << "\n"; + } + } + + Plus_U_Base dftu; + const std::vector l_channel = {2}; + const std::vector hubbard_u = {4.0}; + dftu.init_base(ucell, + 1, // npol + 2, // nspin + l_channel, + false, // yukawa_potential + 0.5, // yukawa_lambda + dir, // global_readin_dir + "", // global_out_dir + "none", // init_chg + "cpu", // device + hubbard_u, + 0.0, // uramping + 1, // init_occ_mat + 0, // mixing_dftu + DFTU_BASE::SOC_LAYOUT_PAULI); + + ASSERT_TRUE(dftu.is_occmat_ready()); + EXPECT_NEAR(dftu.occmat().get(0, 2, 0, 0, 0), occ_up, 1e-12); + EXPECT_NEAR(dftu.occmat().get(0, 2, 1, 4, 4), occ_dn, 1e-12); + + // Before building the potential, uterm_mat is all zeros. + const std::vector>& uterm0 = dftu.get_uterm_mat(); + for (size_t i = 0; i < uterm0.size(); ++i) + { + EXPECT_EQ(uterm0[i].real(), 0.0); + } + + // Build pot_onsite and E_plusU from the file-loaded matrix — this is + // the call the fix adds at the first SCF step for init_occ_mat=1. + pw::compute_pot_uterm_and_energy(ucell, + 2, // nspin + dftu.get_u_current_vec(), + dftu.get_l_channel_vec(), + dftu.get_uterm_mat_index(), + dftu.occmat(), + dftu.get_uterm_mat(), + dftu.energy_ref()); + + const double u_value = 4.0; + const int m_size = 5; + const int size = m_size * m_size; // 25 + const std::vector>& uterm = dftu.get_uterm_mat(); + ASSERT_EQ(uterm.size(), static_cast(2 * size)); + + // spin-up block [0 .. size): diagonal = U*(0.5 - occ_up) = 1.6 + // off-diagonal = U*(0 - occ_up) = -0.4 + for (int m = 0; m < m_size; ++m) + { + EXPECT_NEAR(uterm[m * m_size + m].real(), u_value * (0.5 - occ_up), 1e-12) + << "spin-up diagonal mismatch at m=" << m; + } + EXPECT_NEAR(uterm[0 * m_size + 1].real(), -u_value * occ_up, 1e-12); + EXPECT_NEAR(uterm[1 * m_size + 0].real(), -u_value * occ_up, 1e-12); + + // spin-down block [size .. 2*size): diagonal = U*(0.5 - occ_dn) = 1.2 + // off-diagonal = U*(0 - occ_dn) = -0.8 + for (int m = 0; m < m_size; ++m) + { + EXPECT_NEAR(uterm[size + m * m_size + m].real(), + u_value * (0.5 - occ_dn), 1e-12) + << "spin-down diagonal mismatch at m=" << m; + } + EXPECT_NEAR(uterm[size + 0 * m_size + 1].real(), -u_value * occ_dn, 1e-12); + + // E_U = U * 0.5 * 25 * (occ_up^2 + occ_dn^2) + // = 4.0 * 0.5 * 25 * (0.01 + 0.04) = 2.5 + const double expected_energy = u_value * 0.5 * size + * (occ_up * occ_up + occ_dn * occ_dn); + EXPECT_NEAR(dftu.get_energy(), expected_energy, 1e-12); + EXPECT_GT(dftu.get_energy(), 0.0); +} + +/// For init_occ_mat=2 the occupation matrix is fixed from the file +/// throughout the calculation. Every SCF step must rebuild the potential +/// from the same matrix. This test verifies that (1) the file matrix +/// produces a non-zero potential and energy, and (2) a second call with +/// the unchanged matrix reproduces the same results — the defining +/// property of mode 2 ("fixed matrix"). +TEST_F(DFTUBaseTest, InitOccMat2BuildsPotentialFromFileAndIsStable) +{ + const std::string dir = testing::TempDir(); + const std::string fn = dir + "occ_mat.txt"; + + const double occ_up = 0.1; + const double occ_dn = 0.2; + { + std::ofstream ofs(fn); + ASSERT_TRUE(ofs.is_open()); + ofs << "# compact test fixture\n"; + ofs << " Fe Atom 1 L 2 mag 0.0\n"; + ofs << " spin 1 nelec 0.5\n"; + for (int m0 = 0; m0 < 5; ++m0) + { + for (int m1 = 0; m1 < 5; ++m1) { ofs << " " << occ_up; } + ofs << "\n"; + } + ofs << " spin 2 nelec 0.5\n"; + for (int m0 = 0; m0 < 5; ++m0) + { + for (int m1 = 0; m1 < 5; ++m1) { ofs << " " << occ_dn; } + ofs << "\n"; + } + } + + Plus_U_Base dftu; + const std::vector l_channel = {2}; + const std::vector hubbard_u = {4.0}; + dftu.init_base(ucell, + 1, // npol + 2, // nspin + l_channel, + false, // yukawa_potential + 0.5, // yukawa_lambda + dir, // global_readin_dir + "", // global_out_dir + "none", // init_chg + "cpu", // device + hubbard_u, + 0.0, // uramping + 2, // init_occ_mat — mode 2: fixed matrix + 0, // mixing_dftu + DFTU_BASE::SOC_LAYOUT_PAULI); + + ASSERT_TRUE(dftu.is_occmat_ready()); + + // First call — build potential from the file-loaded matrix. + pw::compute_pot_uterm_and_energy(ucell, + 2, + dftu.get_u_current_vec(), + dftu.get_l_channel_vec(), + dftu.get_uterm_mat_index(), + dftu.occmat(), + dftu.get_uterm_mat(), + dftu.energy_ref()); + + const double energy_first = dftu.get_energy(); + EXPECT_GT(energy_first, 0.0); + + // Snapshot the potential after the first build. + const std::vector> uterm_first = dftu.get_uterm_mat(); + + // Second call — mode 2 means the matrix never changes, so the + // potential and energy must be identical. + pw::compute_pot_uterm_and_energy(ucell, + 2, + dftu.get_u_current_vec(), + dftu.get_l_channel_vec(), + dftu.get_uterm_mat_index(), + dftu.occmat(), + dftu.get_uterm_mat(), + dftu.energy_ref()); + + const double energy_second = dftu.get_energy(); + EXPECT_NEAR(energy_first, energy_second, 1e-12); + + const std::vector>& uterm_second = dftu.get_uterm_mat(); + ASSERT_EQ(uterm_first.size(), uterm_second.size()); + for (size_t i = 0; i < uterm_first.size(); ++i) + { + EXPECT_NEAR(uterm_first[i].real(), uterm_second[i].real(), 1e-12) + << "uterm mismatch at index " << i; + } +} + +/// When DFT+U is off (dft_plus_u <= 0) append_ion_step_snapshot must +/// return immediately, BEFORE the iter >= 1 validation. This guard-first +/// ordering was introduced in step 3/3 of the EXX restart refactor so that +/// a transient iter=0 (e.g. from an upstream EXX rerun signal) does not +/// trigger WARNING_QUIT on a non-DFT+U calculation. Passing iter=0 with +/// dft_plus_u=0 must therefore return cleanly. +TEST_F(DFTUBaseTest, AppendSnapshotGuardDftuOffBeforeIterValidation) +{ + Plus_U_Base dftu; + this->init_dftu(dftu, false); + + // out_occ_mat=true but dft_plus_u=0: the guard must catch dft_plus_u + // and return before the iter < 1 check fires. + const DFTU_BASE::OccmatOutputCfg cfg = {1, 1, 5, true, 0}; + + // Use a unique subdirectory so files from other tests in the same + // binary do not cause false positives. + const std::string out_dir = testing::TempDir() + "guard_test/"; + + // iter=0 is intentionally invalid; the guard-first check ensures + // this does not abort. + DFTU_BASE::append_ion_step_snapshot(dftu, + ucell, + out_dir, + 2, // nspin + 1, // npol + 0, // istep + 0, // iter — invalid, but DFT+U is off + false, + false, + 1e-6, + 0.5, + cfg, + DFTU_BASE::SOC_LAYOUT_PAULI); + + // If we reach here the guard worked. No file should have been written + // (the directory may not even exist, which is fine — the guard + // returned before any file IO). + const std::string fn_g1 = out_dir + "occ_matg1.txt"; + std::ifstream ifs(fn_g1.c_str()); + EXPECT_FALSE(ifs.is_open()); + + // Sanity check: with dft_plus_u > 0 and iter=0 the function WOULD + // abort. We cannot test that path in-process (WARNING_QUIT calls + // exit/MPI_Abort), so the guard-first test above is the best we can do. +} + +// Reading the occupation-matrix file broadcasts on MPI_COMM_WORLD, so +// the test binary must initialize MPI even when ctest launches it as a +// single process. +int main(int argc, char** argv) +{ +#ifdef __MPI + MPI_Init(&argc, &argv); +#endif + testing::InitGoogleTest(&argc, argv); + const int result = RUN_ALL_TESTS(); +#ifdef __MPI + MPI_Finalize(); +#endif + return result; +} diff --git a/tests/01_PW/CASES_CPU.txt b/tests/01_PW/CASES_CPU.txt index bd8cf7e1a63..b0d6cade3e8 100644 --- a/tests/01_PW/CASES_CPU.txt +++ b/tests/01_PW/CASES_CPU.txt @@ -119,8 +119,8 @@ scf_brav_bco_smear_mv 812_PW_LT_sm 813_PW_LT_bacm 814_PW_LT_triclinic -scf_dftu_spin2 -scf_dftu_spin4 +relax_dftu_spin2 +relax_dftu_spin4 scf_deltaspin2 scf_deltaspin4 221_PW_PCC0D_H3O diff --git a/tests/01_PW/CASES_GPU.txt b/tests/01_PW/CASES_GPU.txt index f4dc0654527..fe156d9365b 100644 --- a/tests/01_PW/CASES_GPU.txt +++ b/tests/01_PW/CASES_GPU.txt @@ -111,7 +111,7 @@ scf_brav_bco_smear_mv 812_PW_LT_sm 813_PW_LT_bacm 814_PW_LT_triclinic -scf_dftu_spin2 -scf_dftu_spin4 +#relax_dftu_spin2 +#relax_dftu_spin4 scf_deltaspin2 scf_deltaspin4 diff --git a/tests/01_PW/relax_dftu_spin2/INPUT b/tests/01_PW/relax_dftu_spin2/INPUT new file mode 100644 index 00000000000..d8c31071740 --- /dev/null +++ b/tests/01_PW/relax_dftu_spin2/INPUT @@ -0,0 +1,46 @@ +INPUT_PARAMETERS + +# general +suffix autotest +calculation relax +relax_nmax 3 +basis_type pw +device cpu + +# plane wave +nspin 2 +ecutwfc 10 +nbands 28 +gamma_only 0 +scf_thr 1e-07 +scf_nmax 70 +out_chg 0 +out_occ_mat 1 + +# smearing +smearing_method gaussian +smearing_sigma 0.01 + +# charge mixing +mixing_type broyden +mixing_beta 0.4 +ks_solver dav_subspace +symmetry 0 + +# DFT+U parameters +dft_plus_u 1 +orbital_corr 2 +hubbard_u 5.0 +onsite_radius 3.0 + +init_occ_mat 2 +read_file_dir . +out_freq_elec 1 +out_freq_ion 2 + +pseudo_dir ../../PP_ORB/ +orbital_dir ../../PP_ORB/ +pw_seed 1 + +cal_force 1 +cal_stress 1 diff --git a/tests/01_PW/scf_dftu_spin4/KPT b/tests/01_PW/relax_dftu_spin2/KPT similarity index 68% rename from tests/01_PW/scf_dftu_spin4/KPT rename to tests/01_PW/relax_dftu_spin2/KPT index 35597cecff1..853daf6d9cb 100644 --- a/tests/01_PW/scf_dftu_spin4/KPT +++ b/tests/01_PW/relax_dftu_spin2/KPT @@ -1,4 +1,4 @@ K_POINTS 0 Monkhorst-Pack -2 2 2 0 0 0 +2 1 1 0 0 0 diff --git a/tests/01_PW/scf_dftu_spin2/README b/tests/01_PW/relax_dftu_spin2/README similarity index 100% rename from tests/01_PW/scf_dftu_spin2/README rename to tests/01_PW/relax_dftu_spin2/README diff --git a/tests/01_PW/scf_dftu_spin2/STRU b/tests/01_PW/relax_dftu_spin2/STRU similarity index 100% rename from tests/01_PW/scf_dftu_spin2/STRU rename to tests/01_PW/relax_dftu_spin2/STRU diff --git a/tests/01_PW/relax_dftu_spin2/occ_mat.txt b/tests/01_PW/relax_dftu_spin2/occ_mat.txt new file mode 100644 index 00000000000..51b4a3885a8 --- /dev/null +++ b/tests/01_PW/relax_dftu_spin2/occ_mat.txt @@ -0,0 +1,34 @@ +# ABACUS version: v3.11.0-beta10 +# Written at 2026-10-06 16:07:56 +# RELAX STEP 1 +# Electronic step 25 +# scf_thr 1.00000000e-07 +# drho 9.07140846e-08 + + Fe Atom 1 L 2 mag 1.09213471 + spin 1 nelec 1.81371845 + 0.35878875 -0.03675006 -0.03574056 0.09921064 0.03464868 + -0.03675006 0.43801117 -0.01154315 0.06122753 0.06242945 + -0.03574056 -0.01154315 0.33467807 0.06190739 0.01154211 + 0.09921064 0.06122753 0.06190739 0.24423067 -0.06244040 + 0.03464868 0.06242945 0.01154211 -0.06244040 0.43800979 + spin 2 nelec 0.72158374 + 0.12235312 -0.03335235 0.02008706 -0.05893818 0.04539065 + -0.03335235 0.13260557 -0.09430838 0.07166331 -0.11610061 + 0.02008706 -0.09430838 0.14366414 -0.03479970 0.09429968 + -0.05893818 0.07166331 -0.03479970 0.19037294 -0.06470567 + 0.04539065 -0.11610061 0.09429968 -0.06470567 0.13258798 + + Fe Atom 2 L 2 mag -1.09383536 + spin 1 nelec 0.72078886 + 0.12124536 -0.03351276 0.01918194 -0.05937781 0.04613034 + -0.03351276 0.13491023 -0.09262239 0.07261522 -0.11839744 + 0.01918194 -0.09262239 0.13992032 -0.03322173 0.09262159 + -0.05937781 0.07261522 -0.03322173 0.18980368 -0.06532927 + 0.04613034 -0.11839744 0.09262159 -0.06532927 0.13490927 + spin 2 nelec 1.81462422 + 0.35794004 -0.03700722 -0.03562436 0.09950692 0.03423271 + -0.03700722 0.43859537 -0.01007002 0.06089258 0.06164872 + -0.03562436 -0.01007002 0.33639679 0.06167859 0.01006413 + 0.09950692 0.06089258 0.06167859 0.24309573 -0.06248398 + 0.03423271 0.06164872 0.01006413 -0.06248398 0.43859628 diff --git a/tests/01_PW/relax_dftu_spin2/result.ref b/tests/01_PW/relax_dftu_spin2/result.ref new file mode 100644 index 00000000000..121603b73b8 --- /dev/null +++ b/tests/01_PW/relax_dftu_spin2/result.ref @@ -0,0 +1,5 @@ +etotref -5654.4186834314468797 +etotperatomref -2827.2093417157 +totalforceref 40.305058 +totalstressref 103325.869757 +totaltimeref 2.63 diff --git a/tests/01_PW/relax_dftu_spin4/INPUT b/tests/01_PW/relax_dftu_spin4/INPUT new file mode 100644 index 00000000000..84a3036192d --- /dev/null +++ b/tests/01_PW/relax_dftu_spin4/INPUT @@ -0,0 +1,47 @@ +INPUT_PARAMETERS +suffix autotest +calculation relax +relax_nmax 2 +basis_type pw +device cpu + +# plane wave +noncolin 1 +nspin 4 +ecutwfc 10 +gamma_only 0 +scf_thr 1.0e-7 +scf_nmax 30 +out_chg 0 +out_occ_mat 1 + +# smearing +smearing_method gaussian +smearing_sigma 0.01 + +# charge mixing +mixing_type broyden +mixing_beta 0.4 +ks_solver dav_subspace +symmetry 0 + +# DFT+U parameters +dft_plus_u 1 +orbital_corr 2 +hubbard_u 5.0 +onsite_radius 3.0 + +init_occ_mat 2 +read_file_dir . + +out_freq_elec 1 +out_freq_ion 2 + +kpar 1 +pseudo_dir ../../PP_ORB/ +orbital_dir ../../PP_ORB/ +pw_seed 1 + +out_stru 1 +cal_force 1 +cal_stress 1 diff --git a/tests/01_PW/scf_dftu_spin2/KPT b/tests/01_PW/relax_dftu_spin4/KPT similarity index 68% rename from tests/01_PW/scf_dftu_spin2/KPT rename to tests/01_PW/relax_dftu_spin4/KPT index 35597cecff1..54254473ce8 100644 --- a/tests/01_PW/scf_dftu_spin2/KPT +++ b/tests/01_PW/relax_dftu_spin4/KPT @@ -1,4 +1,4 @@ K_POINTS 0 Monkhorst-Pack -2 2 2 0 0 0 +1 2 1 0 0 0 diff --git a/tests/01_PW/scf_dftu_spin4/README b/tests/01_PW/relax_dftu_spin4/README similarity index 100% rename from tests/01_PW/scf_dftu_spin4/README rename to tests/01_PW/relax_dftu_spin4/README diff --git a/tests/01_PW/scf_dftu_spin4/STRU b/tests/01_PW/relax_dftu_spin4/STRU similarity index 100% rename from tests/01_PW/scf_dftu_spin4/STRU rename to tests/01_PW/relax_dftu_spin4/STRU diff --git a/tests/01_PW/relax_dftu_spin4/occ_mat.txt b/tests/01_PW/relax_dftu_spin4/occ_mat.txt new file mode 100644 index 00000000000..02690cc5f07 --- /dev/null +++ b/tests/01_PW/relax_dftu_spin4/occ_mat.txt @@ -0,0 +1,58 @@ +# ABACUS version: v3.11.0-beta10 +# Written at 2026-10-08 15:33:51 +# RELAX STEP 1 +# Electronic step 30 +# scf_thr 1.00000000e-07 +# drho 5.21682098e-07 + + Fe Atom 1 L 2 mag 1.09308347 0.00423053 -0.01151831 + spin 1 nelec 1.26182180 + 0.23925726 -0.00812736 -0.03495918 -0.02066839 0.03999798 + -0.00812736 0.23817903 -0.05298870 -0.01399639 0.05299584 + -0.03495918 -0.05298870 0.28427727 -0.06639206 -0.02814308 + -0.02066839 -0.01399639 -0.06639206 0.21579949 0.06348580 + 0.03999798 0.05299584 -0.02814308 0.06348580 0.28430876 + spin 12 re + 0.11837214 -0.02789799 -0.00175819 -0.07886795 -0.00530300 + -0.02789799 0.09552847 0.04136487 -0.04828102 -0.04139152 + -0.00175819 0.04136487 0.15267561 0.00517847 0.08931465 + -0.07886795 -0.04828102 0.00517847 0.02729998 -0.00114526 + -0.00530300 -0.04139152 0.08931465 -0.00114526 0.15266553 + spin 12 im + 0.00050911 0.00025731 -0.00017638 0.00106774 0.00036545 + 0.00025731 -0.00050248 -0.00034517 0.00044548 0.00034506 + -0.00017638 -0.00034517 0.00018188 -0.00052221 0.00072431 + 0.00106774 0.00044548 -0.00052221 0.00174285 0.00041333 + 0.00036545 0.00034506 0.00072431 0.00041333 0.00018391 + spin 2 nelec 1.27334011 + 0.24158252 -0.00760473 -0.03513003 -0.01995942 0.03996430 + -0.00760473 0.24006364 -0.05291122 -0.01309021 0.05291653 + -0.03513003 -0.05291122 0.28636121 -0.06644837 -0.02563484 + -0.01995942 -0.01309021 -0.06644837 0.21893609 0.06366294 + 0.03996430 0.05291653 -0.02563484 0.06366294 0.28639665 + + Fe Atom 2 L 2 mag -1.09435018 0.00377374 0.00185987 + spin 1 nelec 1.26873236 + 0.24011019 -0.00795801 -0.03502841 -0.02029743 0.04002410 + -0.00795801 0.23760412 -0.05139792 -0.01386379 0.05140273 + -0.03502841 -0.05139792 0.28729028 -0.06648009 -0.02795862 + -0.02029743 -0.01386379 -0.06648009 0.21644104 0.06361270 + 0.04002410 0.05140273 -0.02795862 0.06361270 0.28728672 + spin 12 re + -0.11843141 0.02731586 0.00171396 0.07943021 0.00600621 + 0.02731586 -0.09834545 -0.04131415 0.04739689 0.04130889 + 0.00171396 -0.04131415 -0.15188266 -0.00592538 -0.08999641 + 0.07943021 0.04739689 -0.00592538 -0.02663553 0.00147915 + 0.00600621 0.04130889 -0.08999641 0.00147915 -0.15188004 + spin 12 im + 0.00027673 -0.00039669 -0.00014033 0.00003725 0.00024354 + -0.00039669 0.00078439 0.00021167 -0.00068912 -0.00021172 + -0.00014033 0.00021167 0.00025385 -0.00036333 0.00018633 + 0.00003725 -0.00068912 -0.00036333 0.00031914 0.00030362 + 0.00024354 -0.00021172 0.00018633 0.00030362 0.00025276 + spin 2 nelec 1.26687249 + 0.23922877 -0.00830010 -0.03543993 -0.01971854 0.04023441 + -0.00830010 0.23889998 -0.05120133 -0.01445721 0.05120779 + -0.03543993 -0.05120133 0.28625648 -0.06696063 -0.02882124 + -0.01971854 -0.01445721 -0.06696063 0.21623666 0.06420774 + 0.04023441 0.05120779 -0.02882124 0.06420774 0.28625058 diff --git a/tests/01_PW/relax_dftu_spin4/result.ref b/tests/01_PW/relax_dftu_spin4/result.ref new file mode 100644 index 00000000000..0d80afa8bc1 --- /dev/null +++ b/tests/01_PW/relax_dftu_spin4/result.ref @@ -0,0 +1,5 @@ +etotref -5652.1157954576747215 +etotperatomref -2826.0578977288 +totalforceref 18.944306 +totalstressref 102121.239260 +totaltimeref 4.24 diff --git a/tests/01_PW/scf_dftu_spin2/INPUT b/tests/01_PW/scf_dftu_spin2/INPUT deleted file mode 100644 index 88bcde220e8..00000000000 --- a/tests/01_PW/scf_dftu_spin2/INPUT +++ /dev/null @@ -1,29 +0,0 @@ -INPUT_PARAMETERS -suffix autotest -calculation scf -basis_type pw -ecutwfc 50 -gamma_only 0 -device cpu - -nspin 2 -nbands 28 -scf_thr 1.0e-6 -scf_nmax 100 -out_chg 0 -smearing_method gaussian -smearing_sigma 0.01 -mixing_type broyden -mixing_beta 0.4 -ks_solver dav_subspace -symmetry 0 - -# DFT+U parameters -dft_plus_u 1 -orbital_corr 2 -hubbard_u 5.0 -onsite_radius 3.0 - -pseudo_dir ../../PP_ORB -orbital_dir ../../PP_ORB -pw_seed 1 diff --git a/tests/01_PW/scf_dftu_spin2/result.ref b/tests/01_PW/scf_dftu_spin2/result.ref deleted file mode 100644 index f9ddfdd28af..00000000000 --- a/tests/01_PW/scf_dftu_spin2/result.ref +++ /dev/null @@ -1,3 +0,0 @@ -etotref -6792.3335167095001452 -etotperatomref -3396.1667583548 -totaltimeref 21.07 diff --git a/tests/01_PW/scf_dftu_spin4/INPUT b/tests/01_PW/scf_dftu_spin4/INPUT deleted file mode 100644 index eb524f797e1..00000000000 --- a/tests/01_PW/scf_dftu_spin4/INPUT +++ /dev/null @@ -1,30 +0,0 @@ -INPUT_PARAMETERS -suffix autotest -calculation scf -basis_type pw -ecutwfc 20 -gamma_only 0 -device cpu - -noncolin 1 -nspin 4 -scf_thr 1.0e-6 -scf_nmax 50 -out_chg 0 -smearing_method gaussian -smearing_sigma 0.01 -mixing_type broyden -mixing_beta 0.4 -ks_solver dav_subspace -symmetry 0 - -# DFT+U parameters -dft_plus_u 1 -orbital_corr 2 -hubbard_u 5.0 -onsite_radius 3.0 - -kpar 1 -pseudo_dir ../../PP_ORB -orbital_dir ../../PP_ORB -pw_seed 1 diff --git a/tests/01_PW/scf_dftu_spin4/result.ref b/tests/01_PW/scf_dftu_spin4/result.ref deleted file mode 100644 index 8242af7627b..00000000000 --- a/tests/01_PW/scf_dftu_spin4/result.ref +++ /dev/null @@ -1,3 +0,0 @@ -etotref -6364.2658763901727070 -etotperatomref -3182.1329381951 -totaltimeref 7.82 diff --git a/tests/17_DS_DFTU/README.md b/tests/17_DS_DFTU/README.md index a3b87cd3d01..02b7d309526 100644 --- a/tests/17_DS_DFTU/README.md +++ b/tests/17_DS_DFTU/README.md @@ -175,9 +175,9 @@ The following test cases are disabled in `CASES_CPU.txt` (commented out with `#` - 09 (PW DFT+U + noncollinear): Only supports **2-process MPI** execution, `result.ref` reference files provided - The following test cases set `kpar=2` in INPUT and require at least **2 MPI processes** to run: 11, 12, 14, 15, 16, 18, 19, 21, 37, 39, 41, 43, 45 -- 62 (LCAO_DFTU_NSCF_Band_XY): Single-thread and multi-thread results are inconsistent; investigation shows HR, HK, and SK are consistent across threads, but eigenvalues from genelpa differ; switching to scalapack_gvx produces consistent results across thread counts. Note: this test is named "NSCF" but actually runs with `calculation = scf` (`scf_nmax = 1`), using pre-shipped charge density and dm_onsite.txt files as initial guess +- 62 (LCAO_DFTU_NSCF_Band_XY): Single-thread and multi-thread results are inconsistent; investigation shows HR, HK, and SK are consistent across threads, but eigenvalues from genelpa differ; switching to scalapack_gvx produces consistent results across thread counts. Note: this test is named "NSCF" but actually runs with `calculation = scf` (`scf_nmax = 1`), using pre-shipped charge density and occ_mat.txt files as initial guess - All NSCF tests (55, 60, 61, 62, 63, 64) have been **converted to SCF+NSCF workflow**: - - Pre-converged `autotest-CHARGE-DENSITY.restart` and `dm_onsite.txt` files have been removed + - Pre-converged `autotest-CHARGE-DENSITY.restart` and `occ_mat.txt` files have been removed - Each test directory contains a `scf/` subdirectory with SCF input files - Run with: `bash ../run_scf_nscf.sh [mpi_np]` - These tests are **disabled in CI** (commented out in CASES_CPU.txt) diff --git a/tests/17_DS_DFTU/run_scf_nscf.sh b/tests/17_DS_DFTU/run_scf_nscf.sh index 7105cdccc4c..35bde1eec2a 100755 --- a/tests/17_DS_DFTU/run_scf_nscf.sh +++ b/tests/17_DS_DFTU/run_scf_nscf.sh @@ -7,7 +7,7 @@ # # This script: # 1. Runs SCF calculation in scf/ subdirectory -# 2. Copies charge density (and dm_onsite.txt for DFT+U) from SCF output +# 2. Copies charge density (and occ_mat.txt for DFT+U) from SCF output # 3. Runs NSCF calculation in current directory # 4. Compares output with reference # @@ -56,11 +56,11 @@ if [ ! -f "${TEST_DIR}/STRU" ]; then exit 1 fi -# Check if this is a DFT+U test (needs dm_onsite.txt handling) +# Check if this is a DFT+U test (needs occ_mat.txt handling) IS_DFTU=false if grep -q "dft_plus_u" "${TEST_DIR}/INPUT" 2>/dev/null; then IS_DFTU=true - echo "DFT+U test detected: will handle dm_onsite.txt files" + echo "DFT+U test detected: will handle occ_mat.txt files" fi # ------------------------------------------------------- @@ -135,15 +135,15 @@ NSCF_CHG_FILE="${NSCF_SUFFIX}-${CHG_BASENAME#*-}" cp "${CHG_FILE}" "${TEST_DIR}/${NSCF_CHG_FILE}" echo " Copied to: ${TEST_DIR}/${NSCF_CHG_FILE}" -# For DFT+U tests, also copy dm_onsite.txt if it exists. -# dm_onsite.txt is read from read_file_dir (global_readin_dir), same as charge density +# For DFT+U tests, also copy occ_mat.txt if it exists. +# occ_mat.txt is read from read_file_dir (global_readin_dir), same as charge density if [ "${IS_DFTU}" = true ]; then - ONSITE_FILE=$(find "${SCF_OUT}" -name "dm_onsite.txt" 2>/dev/null | head -1) + ONSITE_FILE=$(find "${SCF_OUT}" -name "occ_mat.txt" 2>/dev/null | head -1) if [ -n "${ONSITE_FILE}" ]; then - cp "${ONSITE_FILE}" "${TEST_DIR}/dm_onsite.txt" - echo " Copied dm_onsite.txt to: ${TEST_DIR}/dm_onsite.txt" + cp "${ONSITE_FILE}" "${TEST_DIR}/occ_mat.txt" + echo " Copied occ_mat.txt to: ${TEST_DIR}/occ_mat.txt" else - echo " WARNING: dm_onsite.txt not found in SCF output" + echo " WARNING: occ_mat.txt not found in SCF output" fi fi