Ouro DFT now holds gadolinium's 4f shell: what was wrong, and how we built a basis that fixes it

Ouro DFT (ABACUS) could run rare-earth structures, and the runs converged. They were also wrong. This post covers how we found that, why no off-the-shelf fix existed, what it took to build one, and where rare-earth DFT on Ouro stands now. Gadolinium is done. Neodymium and samarium are in progress.

The symptom

While following up the gadolinium literature check, we ran Magnetic moments on hcp Gd (mp-155). It converged normally and reported 0.055 μB per Gd. The measured value is 7.55–7.63 μB. The Mulliken breakdown showed why: −0.05 electrons in the 4f channel and 9.48 in 5d. Gd has seven 4f electrons and about one 5d electron.

The cause: a pseudopotential and a basis that don't belong together

ABACUS runs in a localized basis of numerical atomic orbitals, and each orbital set is generated for one specific pseudopotential. The service paired:

• Pseudopotential: the jingslaw PseudoDojo fully relativistic lanthanide set, with 4f in valence (18 electrons for Gd).

• Orbitals: ABACUS's Dojo-NC-SR_La-Series. We checked the upstream generation inputs: they were built from a scalar-relativistic pseudopotential with 4f frozen in the core. Their one f function is a polarization shell, not a 4f orbital.

A basis with no 4f orbital cannot hold 4f electrons, so SCF pushed them into 5d and the moment collapsed. Nothing crashed and nothing warned. Every lanthanide from Ce to Yb was affected, so the Nd and Sm moments were wrong in the same way. The first thing we shipped was a guard: until a matched basis exists, the service refuses lanthanides with an explanation (example).

Why this isn't a download

The ABACUS ecosystem only distributes lanthanide orbitals with 4f in core: the La-series set above, and the tested APNS bundle, whose lanthanide file is literally named lanthanides-f--core. The 4f-in-core approach is standard and useful, but it has no 4f moment and no 4f anisotropy, which are the two things rare-earth magnets are made of. So we generated our own basis with ABACUS-CSW-NAO. It fits contracted orbitals to reference calculations of Gd dimers at six bond lengths (1.80–5.00 Å) plus the isolated atom, all using the same pseudopotential.

What we learned getting it right

1. Converge the primitive basis before fitting anything. CSW-NAO contracts a large set of primitive spherical waves, and no fit can beat its primitive. We started from the upstream Nd example's 60 Ry primitive cutoff. The first basis held the moment at one volume but had no E(V) minimum and sat 131 eV/atom above plane waves. The isolated atom showed the cause immediately:

Gd's 4f shell is compact, and a coarse primitive can't resolve it.

2. The pseudopotential itself needs 150 Ry. The plane-wave column shows it: the atom moves 13.9 eV from 100 to 150 Ry and 0.06 eV from 150 to 200. The service default of 100 Ry is fine for 3d metals but not here. Lanthanide cells now default to 150 Ry, and explicit lower values are refused.

3. The small basis isn't enough. Checked against plane waves with the same pseudopotential on hcp Gd across a ±6% lattice range:

The plane-wave reference lands on experiment (V₀ ≈ 33.1 ų). The shipped basis carries about 7.4 f electrons with about 7.1 μB of 4f moment. The energy offset to plane waves is 0.18–0.21 eV/atom and drifts about 5 meV/atom rms near V₀, which is normal for an atomic-orbital basis. All numbers are here:

Results on the live service

Same route, same CIF as the failing run (run):

GdCo₅ comes out ferrimagnetic. Heavy rare earths align antiparallel to Fe and Co, and this is exactly where Prophet fails

The order is right. The net moment is lower than the measured ~1.7 μB/f.u., and we know why: these runs are collinear and spin-only. Cobalt carries an orbital moment of about 0.2–0.3 μB per atom in experiment, roughly 1–1.5 μB per formula unit in total. Capturing it needs spin–orbit coupling.

What this means if you run rare-earth DFT on Ouro

• Gd is supported. Other lanthanides are refused until each has a validated basis. Nd and Sm are being generated now.

• Gd cells run at 150 Ry by default, with a basis of 49 orbitals per Gd. Expect them to cost more than an equivalent 3d cell.

• Collinear, plain PBE only for now. No Hubbard U on the 4f shell yet, and 4f magnetocrystalline anisotropy isn't validated. Gd has almost no 4f anisotropy anyway (L = 0), but Nd and Sm are the whole point, and they'll need spin–orbit coupling, DFT+U and control over the 4f occupation. That's the next phase.

• Cached results can't mix bases. Lanthanide cache keys now carry a basis-library revision, so the old 0.06 μB result can never be served again.

Notes for anyone generating bases

• Run the reference geometries in parallel. CSW-NAO runs them one after another. Laying out the folders first, running each geometry in its own container and then fitting cut wall time from 3–5 hours to about 30 minutes, plus about an hour for the fit.

• Pin threads to 1 per MPI rank (OMP_NUM_THREADS, OPENBLAS_NUM_THREADS). Modal sets OMP_NUM_THREADS to the core count at runtime, and 32 ranks × 32 threads ran about 12× slower.

• CSW-NAO needs ABACUS ≥3.7.5 and <3.9.0.6. ABACUS 3.9 renamed OUT.*/kpoints to KPT.info, which CSW-NAO doesn't read. The layout is identical, so a one-line fallback fixes it.

• Keep scratch work on a persistent volume and only resume geometries whose SCF actually converged. Our jobs were preempted five times, and a half-written folder looks finished to a naive resume check.

The generation and validation tools (orbgen_app.py, validate_basis.py) are in the ouro-dft repo. The shipped orbitals sit next to their generation inputs, spillage curves and validation data, so the next element follows the same recipe.