fix(dft+u): stabilize UOM mixing and frozen-occupation potential updates

[liu-687]

Summary

This PR hardens the plane-wave DFT+U occupation-matrix (UOM) feedback loop, with particular focus on noncollinear/SOC calculations using nspin=4. It preserves the existing input interface and user-facing defaults.

For the present SmCo5 system, the source-level changes and the accompanying conservative input settings jointly form the key conditions for stable ordinary-SCF convergence. The source changes make the charge--UOM--potential feedback consistent; the input settings provide a controlled numerical path for introducing the 4f Hubbard interaction. This case-specific convergence contract should not be interpreted as a universal parameter prescription for all 4f materials.

The changes provide:

• an independent, conservative UOM mixer and restart history;
• complete four-component UOM transfer for nspin=4;
• consistent frozen-occupation potential updates for omc=2;
• a guarded Broyden fallback for rank-deficient histories;
• hardened occupation-matrix parsing and MPI-count handling.

The implementation touches only the following source and test files; the numerical validation described below is focused on the PW path:

source/source_base/module_mixing/broyden_mixing.cpp
source/source_base/module_mixing/test/mixing_test.cpp
source/source_estate/module_charge/charge_mixing.cpp
source/source_estate/module_charge/charge_mixing.h
source/source_estate/module_charge/chgmixing.cpp
source/source_estate/test/charge_mixing_test.cpp
source/source_lcao/module_dftu/dftu_lcao_op.cpp
source/source_pw/module_pwdft/dftu_base.cpp
source/source_pw/module_pwdft/dftu_base.h
source/source_pw/module_pwdft/dftu_cal_occ_pw.cpp
source/source_pw/module_pwdft/setup_dftu_pw.cpp

Implementation

Independent UOM mixing

Charge density and the local occupation matrix are different nonlinear variables with different dimensions and numerical stiffness. Reusing the charge-density history for UOM can become unreliable after U-ramping, SCF restart, SOC changes, or occupation-control changes.

The UOM therefore has its own conservative mixer and history. For the reference input, mixing_beta=0.03 gives an effective UOM step of 0.015.

this->mixing_uom = new Base_Mixing::Plain_Mixing(0.5 * this->mixing_beta);

this->mixing_uom->push_data(this->uom_mdata,
                            uom_value_in,
                            uom_value_out,
                            nullptr,
                            true);
this->mixing_uom->cal_coef(this->uom_mdata, inner_product_uom);
this->mixing_uom->mix_data(this->uom_mdata, uom_value_out);

The UOM history is reset together with the charge-mixing history at an SCF restart or U-ramping boundary.

The independent mixer is initialized and reallocated on both the first SCF step and a configured mixing restart:

if (inp.dft_plus_u && inp.mixing_dftu)
{
    p_chgmix->init_mixing_uom();
    dftu.enable_mixing();
    p_chgmix->allocate_mixing_uom(dftu.get_size_pot_uterm_pw());
}

Complete noncollinear UOM transfer

For nspin=4, the UOM contains four Pauli components. The PW path now transfers the complete block instead of only the scalar component:

const int size = (2 * target_l + 1) * (2 * target_l + 1);
const int spin_fold = (this->nspin == 4) ? 4 : 1;

for (int mm = 0; mm < size * spin_fold; ++mm)
{
    this->uom_array[pot_uterm_pw_index[iat] + mm] =
        this->occ_mat[iat][target_l][0][0].c[mm];
}

This keeps longitudinal and transverse spin-orbit components synchronized with the occupation matrix used to construct the DFT+U potential.

Frozen-occupation potential update

omc=2 freezes the occupation matrix, but the DFT+U effective potential must still be rebuilt from that matrix on every SCF step:

if (dftu.get_occ_mat_ctrl() == 2)
{
    dftu.update_eff_pot_pw(ucell);
}
else if (!(iter == 1 && istep == 0))
{
    dftu.cal_occ_pw(psi, wg, ucell, p_chgmix, isk);
}

The ordinary occupation-update path retains its existing first-step behavior. The validation and numerical evidence in this PR are PW-only. A one-line guard is also included in the LCAO occupation lifecycle because the same omc=2 control intentionally freezes the shared occupation matrix:

if ((this->current_spin == this->nspin - 1 || this->nspin == 4)
    && this->dftu->get_occ_mat_ctrl() != 2)
{
    this->dftu->mark_occ_mat_dirty();
}

This guard prevents an intentionally frozen matrix from being marked dirty. It does not change the PW algorithm and is not presented as LCAO numerical validation.

Broyden fallback and input safety

When the small Broyden system is rank deficient, the code now uses the latest linear step, clears the secant history, and emits an explicit warning:

if (info != 0)
{
    ModuleBase::WARNING(
        "Broyden_Mixing",
        "DSYSV failed in Broyden mixing; using the latest linear step.");
    std::fill(coef.begin(), coef.end(), 0.0);
    coef[mdata.start] = 1.0;
    ndim_cal_dF = 0;
    start_dF = -1;
}

Occupation-matrix tags are parsed with full-consumption, range, and atom-bound checks. UOM MPI broadcast sizes are checked before conversion to the int count required by the MPI interface.

std::size_t consumed = 0;
const long parsed = std::stol(number, &consumed, 10);
if (consumed != number.size()
    || parsed < static_cast<long>(std::numeric_limits<int>::min())
    || parsed > static_cast<long>(std::numeric_limits<int>::max()))
{
    return false;
}

if (uom_in.size() > static_cast<std::size_t>(
        std::numeric_limits<int>::max()))
{
    ModuleBase::WARNING_QUIT(
        "Charge_Mixing::mix_uom",
        "UOM buffer is too large for MPI broadcast count");
}
Parallel_Common::bcast_double(
    uom_in.data(), static_cast<int>(uom_in.size()));

Here, “range” refers to the representable integer range and the explicit atom-index bound; the PR does not claim semantic validation of every possible L, ORBITAL, or spin value. The Broyden fallback is also not a convergence declaration: subsequent SCF residual, energy, and scf_nmax criteria remain authoritative.

Stage1 validation case

The reproducible validation case is a six-atom SmCo5 cell with Sm 4f states in the valence, ordinary PW SCF, DFT+U, noncollinear magnetism, and SOC. The relevant settings are:

basis_type       pw
scf_thr_type     1
scf_thr          1e-4
scf_ene_thr      1e-3
scf_nmax         800
mixing_type      plain
mixing_beta      0.03
mixing_beta_mag  0.03
mixing_gg0       0
mixing_gg0_mag   0
mixing_restart   3e-3
mixing_dftu      1
uramping         0.25
omc              0
sc_mag_switch    0
noncolin         1
nspin            4
lspinorb         1
onsite_radius    3.0
dft_plus_u       1
orbital_corr     -1 3
hubbard_u        0 4.7
kpar             4
ecutwfc          220
kspacing         0.14

Full Stage continuation

The case uses a four-stage continuation so that the localized 4f basin is formed first, then tested under fixed U, and finally released to an ordinary SCF. The state transfer is limited to the physical SCF state and occupation matrix:

*-CHARGE-DENSITY.restart  -> next stage charge state through init_chg=file
dm_onsite.txt             -> next-stage occupation-matrix input

Stage	Key control	Role and result
Stage1	omc=0, uramping=0.25, U: 0 -> 4.7 eV	Ordinary PW SCF forms the localized Sm 4f basin; the original reference run used for physical comparison required 351 steps and 20 restart/U-update events.
Stage2	omc=2, uramping=0, fixed U=4.7 eV	Freezes the converged occupation matrix and tests the fixed-U potential. The matrix is frozen, but V_U[n_frozen] is rebuilt at every SCF step.
Stage3	omc=1, uramping=0, fixed U=4.7 eV	Releases the occupation constraint while retaining the fixed-U basin.
Stage4	omc=0, uramping=0, fixed U=4.7 eV	Returns to ordinary PW SCF; the final fixed-U calculation reaches #SCF IS CONVERGED# in 28 steps.

This continuation demonstrates the purpose of the frozen-occupation potential update: omc=2 freezes the occupation matrix, not the corresponding DFT+U potential. Rebuilding V_U from the frozen matrix keeps the Hamiltonian consistent and makes Stage2 a valid fixed-U bridge between basin formation and the final ordinary SCF.

Same-input before/after comparison

The parent develop executable and the PR head were run with the identical input and the same 4-node/128-CPU configuration. The table reports total-cell magnetization, not a direct Sm local-moment measurement; the local Sm 4f UOM signature is given separately below.

Case	Steps	Restart/U events	Final DRHO	Final energy (eV)	Total-cell M_z	SCF status
Parent develop	487	20/20	9.9804e-5	-23259.936324	1.99297912	converged
PR head	347	20/20	3.4190e-5	-23265.502500	1.87015681	converged

The PR head uses 140 fewer electronic steps and reduces the elapsed wall time from approximately 20 h 09 min to 15 h 33 min for these matched runs. Both runs reach U=4.7 eV through the same restart/U-ramping schedule. The PR energy agrees with the previously successful Stage1 reference (-23265.502096 eV) to approximately 4.0e-4 eV; the approximately 5.5662 eV/cell difference from the parent shows that the two runs reached different self-consistent states. This supports improved convergence behavior and consistency with the known reference basin, but the energy difference alone is not a proof of absolute physical superiority.

The 351-step value above is from the original successful Stage1 reference run. The 347-step PR row is a separate matched rerun with the same logical input, so these values should not be read as contradictory deterministic step counts.

The comparison shows why the source changes and input settings must be considered together for this case: independent UOM history and complete nspin=4 state transfer protect the self-consistent map, while plain mixing, conservative mixing strengths, mixing_restart=3e-3, and uramping=0.25 keep the iteration within a controllable basin. The settings are numerical controls around a fixed physical model, not substitutes for the source-level corrections.

4f occupation and magnetic signature

For Sm atom 6 (l=3, ORBITAL=0), the validated Stage1 UOM reference has the following Pauli-component traces:

T0 =  5.39868025
T1 =  0.00000036
T2 = -0.00000010
T3 = -5.36866830

The scalar component eigenvalues are:

0.96019216  0.95966444  0.92161841  0.53848426
0.95575281  0.11064831  0.95231985

Five eigenvalues are close to one, the transverse traces are negligible, and the longitudinal component dominates. This is consistent with a localized Sm 4f5-like basin. The corresponding local UOM magnetic signature is:

Magnetism for atom 6: 0.00000036 -0.00000010 -5.36866830

The component signs follow the implementation's Pauli convention; T3<0 does not represent a negative occupation. The local spin density is primarily aligned with the z axis. These values are a validated UOM signature, not a direct before/after table of Sm local spin and orbital moments; the quantitative total magnetic moment additionally requires the orbital contribution.

Focused validation

• The independent-UOM and full charge-mixing tests pass.
• The singular-history Broyden test verifies the latest-step fallback and explicit DSYSV warning.
• In the matched nspin=4, omc=2 control, the parent reaches the 40-step limit with DRHO=3.4010e-4, while the PR run reaches #SCF IS CONVERGED# at step 39 with DRHO=9.5907e-5.
• An ordinary collinear DFT+U case with nspin=2, no SOC, and no U-ramping converges successfully with the PR executable. This is a compatibility check, not a quantitative develop-versus-PR performance comparison.
• Compact occupation-matrix headers and malformed or out-of-range atom tags are handled with explicit validation diagnostics.

No new user-facing parameter is introduced and no existing default is changed. The PR makes the UOM feedback state explicit and robust while preserving the established DFT+U input interface. For the specified SmCo5 U, onsite_radius, SOC, and magnetic initialization, the matched Stage1 and continuation results support a stable convergence path and consistency with the validated 4f basin. They do not claim universal parameter optimality, a complete physical DSYSV-recovery demonstration, or a quantitative collinear performance guarantee for every 4f system.

[Copilot]

Pull request overview

This PR improves robustness of the plane-wave DFT+U SCF path by decoupling occupation-matrix mixing from charge mixing, ensuring frozen occupation-matrix control (omc=2) still supplies a valid DFT+U effective potential, making nspin=4 occupation-matrix packing consistent, accepting compact occupation-matrix headers, and preventing SCF aborts when Broyden’s secant history becomes singular.

Changes:

• Add an independent, conservative Plain mixer + independent history for the DFT+U occupation matrix, and rebuild it on SCF mixing restart.
• For omc=2, rebuild the PW DFT+U effective potential from the frozen occupation matrix each SCF step; keep the historical “skip first-step occ update” behavior for the normal path.
• Improve robustness/consistency: pack all 4 Pauli blocks for nspin=4, accept Atom=6/L=3/ORBITAL=0 compact headers, and fall back gracefully when DSYSV fails in Broyden mixing (plus add focused unit tests).

Reviewed changes

Copilot reviewed 11 out of 11 changed files in this pull request and generated 2 comments.

Show a summary per file

File	Description
source/source_pw/module_pwdft/setup_dftu_pw.cpp	Rebuild PW DFT+U effective potential for frozen occupation matrices (omc=2) before Hamiltonian construction.
source/source_pw/module_pwdft/dftu_cal_occ_pw.cpp	Add update_eff_pot_pw() helper and ensure nspin=4 UOM packing/sync covers all 4 Pauli blocks.
source/source_pw/module_pwdft/dftu_base.h	Declare the new update_eff_pot_pw() interface for PW DFT+U.
source/source_pw/module_pwdft/dftu_base.cpp	Accept both spaced and compact tagged headers while reading occupation-matrix files.
source/source_lcao/module_dftu/dftu_lcao_op.cpp	Avoid marking the LCAO occupation matrix dirty when omc=2 intentionally freezes it.
source/source_estate/test/charge_mixing_test.cpp	Add a unit test asserting DFT+U occupation-matrix mixing uses an independent mixer and conservative step.
source/source_estate/module_charge/chgmixing.cpp	Initialize the UOM mixer and allocate independent UOM mixing data on first iter and mixing restart.
source/source_estate/module_charge/charge_mixing.h	Add init_mixing_uom() API and store an independent mixing_uom pointer.
source/source_estate/module_charge/charge_mixing.cpp	Implement independent UOM mixer lifecycle and mixing path (separate history/inner product + MPI sync).
source/source_base/module_mixing/test/mixing_test.cpp	Add a unit test for Broyden fallback behavior when the secant history is singular.
source/source_base/module_mixing/broyden_mixing.cpp	Replace abort-on-DSYSV-failure with a “latest linear step” fallback and history reset.

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[mohanchen]

Thanks for the detailed write-up. The source-level changes look reasonable, but before we can make a decision on merging, I'd like to see a systematic before/after comparison rather than only the post-fix result.

Right now the PR presents the Stage1 case as "351 steps, converged, DRHO = 3.6e-5, E = -23265.5 eV", but there is no baseline from the current develop branch on the same input. Without that baseline, we can't actually tell whether the changes are an improvement, neutral, or even a regression:

Convergence behavior: Does develop also converge on this SmCo5 + SOC + DFT+U + U-ramping case? If yes, in how many steps, and to what residual/energy? If the patched version takes 351 steps but develop takes 250, that's a regression in convergence speed, not a stability win.
Divergence/crash cases: If the motivation is "stabilize UOM mixing", there should be at least one concrete input where develop diverges or aborts (e.g. DSYSV non-zero exit, UOM history corruption after restart, nspin=4 inconsistent potential). Pointing at that input and showing the patched version recovers is the actual evidence of stability improvement.
Correctness validation: For the nspin=4 full four-component UOM transfer and the omc=2 potential rebuild, please provide a case where the old behavior produced a physically inconsistent result (e.g. wrong local magnetic moment, non-converged energy) and the new behavior fixes it. A unit test asserting the code path is taken is necessary but not sufficient — it doesn't prove the physics is right.
Separately, even when the intent is a correctness fix, there is real risk of regression or new bugs:

Switching UOM from the shared charge mixer to an independent Plain_Mixing(0.5 * mixing_beta) changes the mixing strength for every DFT+U user, not just the SOC/U-ramping subset. We should verify that ordinary collinear DFT+U cases (nspin=2, no SOC, no U-ramping) do not slow down or lose convergence on existing tests.
The Broyden fallback now silently continues with the latest step and clears the secant history (ndim_cal_dF = 0). This avoids the abort, but if the singularity is a symptom of genuine SCF divergence rather than numerical linear dependence, masking it could let a divergent run drag on to scf_nmax instead of failing fast. Please confirm the warning is still loud enough that users/reviewers can distinguish "recovered" from "papered over".
Could you add a before/after table (steps, DRHO, energy, local magnetic moment, and whether #SCF IS CONVERGED# is reached) for at least the Stage1 input, and ideally one case each exercising the nspin=4, omc=2, and Broyden-singular paths? Once we can see the delta, we can make an informed call on whether this is a net improvement worth merging.

[liu-687]

@mohanchen Thanks for the review. I reran the case on the parent develop executable and on the PR head with the identical Stage1 input and the same 4-node/128-CPU layout.

Case	Electronic steps	Restart/U events	Final DRHO	Final energy (eV)	Total-cell M_z	SCF status
Parent develop	487	20/20	9.9804e-5	-23259.936324	1.99297912	converged
PR head	347	20/20	3.4190e-5	-23265.502500	1.87015681	converged

The PR head uses 140 fewer electronic steps and the elapsed time for these matched runs decreases from about 20 h 09 min to 15 h 33 min. Both runs reach U=4.7 eV through the same 20 restart/U-ramping events, so the difference is not due to a shorter ramp schedule. The approximately 5.5662 eV/cell energy difference shows that the two calculations reached different self-consistent states. The PR energy agrees with the previously successful Stage1 reference (-23265.502096 eV) within approximately 4.0e-4 eV, which supports consistency with that reference basin; the energy difference alone is not treated as proof of absolute physical superiority.

The original successful Stage1 reference required 351 steps. The 347-step value above is from a separate PR-head rerun used for the matched comparison, with the same logical input; the two counts are therefore not expected to be identical.

For the frozen-occupation path, the fixed-U continuation completed through omc=2 -> omc=1 -> omc=0, with the final ordinary omc=0 stage converging in 28 steps. This tests that omc=2 freezes the occupation matrix while rebuilding the corresponding DFT+U potential at every SCF step. It is a convergence/control-path result rather than a standalone proof of every physical observable.

Additional controlled results:

• In the matched nspin=4, omc=2 control, the parent reached the 40-step limit with DRHO=3.4010e-4, while the PR run reached #SCF IS CONVERGED# at step 39 with DRHO=9.5907e-5.
• An ordinary collinear DFT+U case with nspin=2, no SOC, and no U-ramping converged in 28 steps with the PR executable. This is a compatibility check; I have not yet made a quantitative develop-versus-PR performance comparison for that case.
• Compact occupation-matrix headers are accepted, and malformed or out-of-range atom tags terminate through explicit diagnostics. The parser change specifically validates complete integer tokens, representable integer range, and atom bounds; it does not claim semantic validation of every possible L, ORBITAL, or spin value.
• The singular-history unit test confirms the latest-step fallback, secant-history reset, and explicit DSYSV failed in Broyden mixing warning. I have not reproduced a complete physical SCF run in which parent develop fails at DSYSV and the PR version subsequently recovers and converges. The fallback therefore should not be read as a convergence guarantee; the usual SCF residual, energy, and scf_nmax criteria remain authoritative.

For the physical reference, the validated Sm atom-6 4f UOM has

T0 =  5.39868025
T1 =  0.00000036
T2 = -0.00000010
T3 = -5.36866830

The scalar-component eigenvalues are

0.96019216  0.95966444  0.92161841  0.53848426
0.95575281  0.11064831  0.95231985

Five eigenvalues are close to one, one is partially occupied, and the transverse traces are negligible. This is consistent with a localized 4f5-like basin whose spin density is primarily aligned with the z axis. These are UOM occupation/signature diagnostics, not a direct before/after table of Sm local spin and orbital moments. The current nspin=4 local-output routine repeats the scalar block when printing components, so repeated component values are not used as independent local-moment vectors.

Overall, the data support improved convergence behavior and consistency with the previously validated Sm 4f basin for this target case. They do not claim a universal parameter prescription, a complete end-to-end physical DSYSV recovery demonstration, or a quantitative collinear performance guarantee without the corresponding develop baseline.