Fix 8084: DFT+U refactor step 12

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.

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

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

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

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_;
 };

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<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
 	const std::vector<std::vector<TK>>& dm_vec = this->dmat.dm->get_dmk_vec();
 
     // 1) calculate the local occupation number matrix and energy correction in DFT+U
-    finish_dftu_lcao<TK>(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<TK>(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<TK, TR>::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)

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;

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,14 +222,24 @@ void ESolver_KS_PW<T, Device>::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_,
                                  this->stp.template get_psi_t<T, Device>(),
                                  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<T, Device>::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<T, Device>(), 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<T, Device>(),
+                            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);
 }

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<double> energies_Ry;
     std::vector<double> 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);
 

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);

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<char>(ifs)), std::istreambuf_iterator<char>());
-    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"));

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.