SkalaUKS does not converge to spin-symmetry-broken solutions; PBE-UKS converges from the same initial guess

[xieziren]

Summary

For stretched H₂, SkalaUKS fails to converge to the spin-symmetry-broken (open-shell singlet) solution, at every geometry past the point where such a solution becomes the physical one. PBE-UKS converges at all of those geometries from the same broken-symmetry initial guess, same grid, same code path. retry_scf runs to exhaustion without helping.

Broken-symmetry UKS is the standard way to get a qualitatively correct dissociation curve out of a density functional, so this closes off a fairly common use case.

Reproduction

repro.sh builds a clean venv and runs repro.py (both attached below) on a single geometry, about 8 minutes end to end:

H2 at R = 3.0 A, cc-pvtz, grids.level = 5
torch 2.13.0

Reference: does a broken-symmetry solution exist at this geometry?
  UHF, plain           E=   -0.831767  <S^2>=0.000  converged=True
  UHF, stability-follow E=   -0.999771  <S^2>=0.995  converged=True

Control: PBE-UKS from that density matrix
  PBE-UKS              E=   -0.999767  <S^2>=0.989  converged=True

Subject: SkalaUKS from the same density matrix, through retry_scf
  SkalaUKS             E=   -0.940861  <S^2>=0.000019  converged=False
  retry_scf attempts   ntries=5  cycles=203
  min HOMO-LUMO gap    up=-0.06048927499339091  down=-0.04483770898210265

A broken-symmetry solution exists here and PBE-UKS reaches it from the given guess; SkalaUKS ends up back on the closed-shell branch with converged=False after the retry ladder is exhausted.

The scan below is the same procedure over seven bond lengths.

H₂ / cc-pVTZ / grids.level = 5 / with_dftd3=False. Every UKS run starts from the density matrix of a converged, stability-followed UHF solution. Skala runs go through skala.pyscf.retry.retry_scf.

R (Å)	UHF ⟨S²⟩	PBE-UKS: E (Ha)	⟨S²⟩	SkalaUKS: E (Ha)	⟨S²⟩
0.74	0.000	−1.166072	0.000 ✅	−1.173358	0.000 ✅
1.20	0.000	−1.108157	0.000 ✅	−1.108190	0.000 ✅
1.60	0.686	−1.044745	0.000 ✅	−1.045171	0.000 ✅
2.00	0.903	−1.011884	0.712 ✅	−1.001167	0.003 ❌
2.40	0.970	−1.002740	0.925 ✅	−0.270637	0.999 ❌
3.00	0.995	−0.999767	0.989 ✅	−0.940861	2e−05 ❌
4.00	1.000	−0.999272	1.000 ✅	−0.925807	1e−04 ❌

✅/❌ is mf.converged. Totals: UHF 7/7, PBE-UKS 7/7, SkalaUKS 3/7.

The 3/7 overstates it. All three Skala successes are ⟨S²⟩ = 0 closed-shell solutions, and at R ≤ 1.6 PBE returns closed-shell too — the Coulson–Fischer point for these functionals sits between 1.6 and 2.0 Å, so closed-shell is simply the right answer there and I would not count it either way.

The number that matters is the bottom four rows: on the geometries where a broken-symmetry solution exists and PBE-UKS finds it, SkalaUKS reaches it 0 out of 4 times.

retry_scf does not help

I ran the official ladder rather than hand-tuning, so this is not a case of not having tried the supported path. On all four failing geometries it ran to the end and still returned converged=False, with SCFState reporting ntries = 5 and 203 cycles in total. Wall time 6–14 minutes per geometry.

One remark on the ladder's design, offered as a question rather than a claim. The level_shift rung only accepts a result if the gap actually opened:

sufficient_gap = all(gap is not None and gap >= SMALL_GAP for gap in gaps)
scf.converged = scf.converged and sufficient_gap

For a near-degenerate broken-symmetry singlet, opening the HOMO–LUMO gap is more or less the definition of collapsing back onto the closed-shell solution, so this heuristic may be pulling in the wrong direction for exactly this class of system. My R = 3.0 and R = 4.0 runs both end at ⟨S²⟩ ≈ 0, i.e. back on the closed-shell branch. Is there a recommended way to handle spin-symmetry-broken cases, or would an ⟨S²⟩ check in the ladder make sense?

Not an initial-guess artifact, and not the sorted-grid problem

Two things worth ruling out up front, since both have come up before.

The custom initial guess is not what breaks it. I saw #76, where a CP2K convergence problem turned out to come from restarting off a previous wavefunction and went away with RESTART=FALSE. That is not what is happening here. Running SkalaUKS from the default guess at the same four geometries, no custom density matrix anywhere:

R (Å)	E (Ha)	⟨S²⟩	converged
2.00	−1.001037	0.000000	True
2.40	−0.965562	0.000001	True
3.00	−0.941665	0.000000	True
4.00	−0.926783	0.000000	False

So at R = 2.0–3.0 the solver is perfectly happy from a default start — it just lands on the closed-shell branch, which is the wrong branch at these geometries. The failure is specific to reaching the broken-symmetry solution, and a broken-symmetry guess is the only way to ask for it. Supplying one is not an exotic thing to do, and PBE-UKS handles the same guess without complaint.

R = 4.00 fails with no custom guess at all, and it still fails after the full ladder: retry_scf returns converged=False with ntries=5, E = −0.925091, ⟨S²⟩ = 1e−04, min gap −0.0100 / −0.0094. That one needs no initial-guess discussion whatsoever.

Not the unsorted-grid issue from #62 either. mf.grids is UnsortableGrids with alignment = 1 after __init__, and stays that way after kernel() — I set grids.level after construction, and initialize_grids handles that case. The failure mode here is also different: #62 surfaces as LinAlgError: Eigenvalues did not converge, whereas these runs complete every cycle and simply report converged=False.

Two secondary observations

Diagnosed gaps go negative, by implausible amounts. From SCFState, the minimum HOMO–LUMO gap over the iteration:

R (Å)	min gap, up	min gap, down
2.00	−2732.92	−115.87
2.40	−1781.87	−967.67
3.00	−0.060	−0.045
4.00	−0.00097	−0.00017

A failed run can look spin-correct while the energy is meaningless. At R = 2.40 the result is E = −0.2706 Ha, against −1.0064 from FCI in the same basis — off by 462 kcal/mol — while ⟨S²⟩ = 0.9986 looks exactly like a healthy open-shell singlet. Anyone filtering on ⟨S²⟩ rather than on converged would accept that point.

A possible lead — unverified

While looking for a cause I noticed SkalaNumInt.gen_response in skala/pyscf/numint.py selects the Coulomb response branch on ks.mol.spin:

vj = ks.get_j(ks.mol, dm1, hermi=1)
if ks.mol.spin == 0:
    v1 += vj                # restricted convention
else:
    v1 += vj[0] + vj[1]     # unrestricted convention

PySCF's own unrestricted response uses v1 += vj[0] + vj[1] (pyscf/scf/_response_functions.py:212). Branching on mol.spin rather than on whether the object is unrestricted means the two disagree in exactly one window: an unrestricted calculation with mol.spin == 0 — which is precisely the broken-symmetry singlet case. The shapes stay compatible, so nothing raises.

gen_response feeds newton() and stability(), and SOSCF is the last rung of retry_scf, so there is a plausible path from this to the behaviour above.

I have not verified that this is the cause — I have only read the code and noticed the mismatch. It could be unrelated. I mention it in case it saves you time; please treat the convergence data above as the actual report.

What I am not claiming

• Not that the Skala functional cannot describe static correlation. The network is explicitly spin-symmetrised, and this report is about SCF solver behaviour, which is a separate question.
• Not anything derived from the non-converged energies. They are listed to show the failure mode, not as results.
• SkalaRKS is unaffected — it converges at all 21 points of this dissociation curve under the same settings, without needing retry_scf. Genuine high-spin open shells are also fine: CH₂ triplet, O₂ triplet and the H atom all converge with correct ⟨S²⟩.

Environment

Reproduced in a fresh venv containing nothing but skala==2026.8 and its dependencies (61 packages), built by repro.sh:

skala (pip):   2026.8          # skala.__version__ reports '1.1.1'
functional:    skala-1.1
pyscf:         2.13.1
torch:         2.13.0
e3nn:          0.6.0
python:        3.11.15
platform:      macOS arm64 (Apple M4 Max), CPU

Nothing is patched or pinned in that environment. I also get bit-identical numbers — energies, ⟨S²⟩, ntries, cycle count, and both gap values to full precision — from a separate environment on e3nn==0.4.4, so the behaviour does not depend on the e3nn version and is deterministic rather than a sensitivity to the initial guess.

As a footnote on the version string: the distribution installs as 2026.8 while skala.__version__ reports 1.1.1. Not important here, but it makes bug reports harder to pin down.

Files

Save these two side by side and run repro.sh. It builds a throwaway venv, asserts nothing in it is patched, and runs repro.py.

repro.sh

#!/usr/bin/env bash
# Clean-environment reproduction: SkalaUKS does not reach the broken-symmetry
# solution for stretched H2, while PBE-UKS does from the same initial guess.
#
# Builds a throwaway venv, checks that nothing in it is patched, and runs
# repro.py. Save this file and repro.py side by side and run this one.
#
# e3nn >= 0.5 is requested deliberately: under e3nn 0.4.4 with torch >= 2.6,
# `import skala.pyscf` fails outright and the usual workaround is a sitecustomize
# or .pth shim calling torch.serialization.add_safe_globals([slice]). Pinning
# e3nn >= 0.5 avoids needing any such shim, so the environment reported in the
# issue is a stock install. The shim check below asserts that -- if you have one
# installed globally it would silently apply here too, and the point of this
# script is that nothing is patched.

set -euo pipefail

HERE="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
VENV="${1:-$HOME/ai4s-envs/skala-uks-repro}"
PY="$VENV/bin/python"

echo "==> creating clean venv at $VENV"
uv venv "$VENV" --python 3.11

echo
echo "==> installing skala 2026.8 with e3nn >= 0.5 (no shim needed)"
uv pip install --python "$PY" "skala==2026.8" "e3nn>=0.5"

echo
echo "==> confirming the environment is stock (no add_safe_globals shim present)"
SITE="$("$PY" -c 'import site; print(site.getsitepackages()[0])')"
if ls "$SITE"/ai4s_e3nn_fix.py "$SITE"/zz_ai4s.pth >/dev/null 2>&1; then
  echo "    FAIL: torch.serialization shim found in $SITE -- not a clean environment"
  exit 1
fi
echo "    OK: no torch.serialization shim in $SITE"

echo
echo "==> import must succeed with no patching"
"$PY" -c "from skala.pyscf import SkalaUKS; import e3nn, torch; print('OK, e3nn', e3nn.__version__, '| torch', torch.__version__)"

echo
echo "==> recording versions"
uv pip freeze --python "$PY" > "$HERE/versions.txt"

echo
echo "==> running the reproduction (about 8 minutes, mostly the SOSCF rung)"
"$PY" -u "$HERE/repro.py"

echo
echo "==> done. $VENV can be deleted; this script rebuilds it."

repro.py

"""Minimal reproducer: SkalaUKS does not reach the broken-symmetry solution.

Stretched H2 has a spin-symmetry-broken UKS solution. This script:
  1. converges UHF and follows its stability analysis to get that solution,
  2. feeds the resulting density matrix to PBE-UKS   -> converges (control),
  3. feeds the same density matrix to SkalaUKS       -> does not converge,
     even after skala's own retry_scf ladder runs to exhaustion.

Runtime: roughly 8 minutes on a laptop CPU, almost all of it in the SOSCF stage
of the retry ladder for the Skala run.

    python repro.py
"""

import numpy as np
import torch
from pyscf import dft, gto, scf

from skala.pyscf import SkalaUKS
from skala.pyscf.retry import retry_scf

R = 3.0            # Angstrom; well past the DFT Coulson-Fischer point
BASIS = "cc-pvtz"
GRID_LEVEL = 5


def broken_symmetry_uhf(mol):
    """Plain UHF collapses onto the RHF solution here; stability following is required."""
    mf = scf.UHF(mol)
    mf.kernel()
    print(f"  UHF, plain           E={mf.e_tot:12.6f}  <S^2>={mf.spin_square()[0]:.3f}  "
          f"converged={mf.converged}")
    mo = mf.stability()[0]
    mf.kernel(mf.make_rdm1(mo, mf.mo_occ))
    print(f"  UHF, stability-follow E={mf.e_tot:12.6f}  <S^2>={mf.spin_square()[0]:.3f}  "
          f"converged={mf.converged}")
    return mf


def main():
    print(f"H2 at R = {R} A, {BASIS}, grids.level = {GRID_LEVEL}")
    print(f"torch {torch.__version__}\n")

    mol = gto.M(atom=f"H 0 0 0; H 0 0 {R}", basis=BASIS, verbose=0)

    print("Reference: does a broken-symmetry solution exist at this geometry?")
    uhf = broken_symmetry_uhf(mol)
    dm0 = uhf.make_rdm1()
    print()

    print("Control: PBE-UKS from that density matrix")
    pbe = dft.UKS(mol)
    pbe.xc = "pbe"
    pbe.grids.level = GRID_LEVEL
    pbe.max_cycle = 100
    pbe.kernel(dm0)
    print(f"  PBE-UKS              E={pbe.e_tot:12.6f}  <S^2>={pbe.spin_square()[0]:.3f}  "
          f"converged={pbe.converged}\n")

    print("Subject: SkalaUKS from the same density matrix, through retry_scf")
    skala = SkalaUKS(mol, xc="skala-1.1", with_dftd3=False)
    skala.grids.level = GRID_LEVEL
    skala.max_cycle = 100
    # retry_scf calls scf.kernel() with no dm0, so the guess is injected here to make
    # sure every rung of the ladder starts from the broken-symmetry density.
    skala.get_init_guess = lambda *a, **k: dm0
    skala, state = retry_scf(skala)
    gap_up, gap_dn = state.get_gap()
    print(f"  SkalaUKS             E={skala.e_tot:12.6f}  <S^2>={skala.spin_square()[0]:.6f}  "
          f"converged={skala.converged}")
    print(f"  retry_scf attempts   ntries={state.ntries}  cycles={len(state.cycles)}")
    print(f"  min HOMO-LUMO gap    up={gap_up}  down={gap_dn}")

    print("\nExpected on skala 2026.8 / pyscf 2.13.1:")
    print("  UHF and PBE-UKS both converge to <S^2> ~ 0.99;")
    print("  SkalaUKS returns converged=False after the ladder is exhausted.")


if __name__ == "__main__":
    main()

[JensWehner]

Thanks for the detailed description, we will look into that.

[esoteric-ephemera]

It seems to be a mild grid sensitivity issue. Copying the density matrix from the previous separation does seem to help a bit as well, that's typically the approach I would take for these kinds of problems, as SCF stability is hard for higher-level functionals:

from tqdm import tqdm

from skala.pyscf.dft import SkalaUKS
from pyscf.gto import M
from skala.functional import load_functional
import numpy as np

def run_h2_dissoc(sep : float, atom : bool = False, dm : np.ndarray | None = None) -> tuple[float,np.ndarray]:
    mol = M(
        atom = f"H 0 0 0; H 0 0 {sep}" if not atom else f"H 0 0 0",
        basis = "cc-pvtz",
        spin = 0 if not atom else 1,
        unit = "Angstrom",
        charge = 0,
    )
    scf = SkalaUKS(mol,xc=load_functional("skala-1.1"))
    scf.max_cycle = 100
    scf.grids.level = 3
    toten = scf.kernel(dm0 = dm)
    if not scf.converged:
        raise RuntimeError("SCF did not converge.")
    return toten, scf.make_rdm1()

atom_energy,_ = run_h2_dissoc(0.,atom=True)
seps = np.linspace(0.5, 4, 10)
dissoc_no_dm = np.zeros(seps.shape)
dissoc_from_dm = np.zeros(seps.shape)
dm = None
for i, sep in tqdm(enumerate(seps),total=len(seps)):
    try:
        toten, _ =  run_h2_dissoc(sep)
        dissoc_no_dm[i] = toten - 2*atom_energy
    except RuntimeError:
        dissoc_no_dm[i] = np.nan

    try:
        toten, dm =  run_h2_dissoc(sep,dm=dm)
        dissoc_from_dm[i] = toten - 2*atom_energy
    except RuntimeError:
        dissoc_from_dm[i] = np.nan

With Grid level 3 in pyscf, here's what I find (nan indicated non-convergence). For separations more than twice the Coulson-Fischer point, there are probably numerical issues with basically treating two dissociated non-overlapping densities:

Separation (Å)	Dissociation Energy (Default init)	Dissociation Energy (Init from previous DM)
0.500000	-0.102332	-0.102332
0.888889	-0.162352	-0.162352
1.277778	-0.094663	-0.094663
1.666667	-0.036790	-0.036790
2.055556	0.003762	0.003762
2.444444	0.037678	0.036991
2.833333	0.053036	0.053036
3.222222	0.063304	nan
3.611111	nan	0.069398
4.000000	nan	nan