Ouro DFT's MAE route: an idle A100, a slow TB2J loop, and a band cut that made results depend on the cell

Magnetic anisotropy energy on Ouro DFT (ABACUS) went through a round of benchmarking this week. That turned up three problems. Two made MAE jobs slow and expensive. The third was a numerical bias in the result itself. All three are fixed on the live service. This post covers what we found, how we checked it, and what changes for anyone using MAE numbers.

Short version:

• MAE jobs now run on CPUs and are 2–4× faster and 6–12× cheaper.

• The FePt MAE drops from 3.16 to 2.87 meV/f.u. It now matches self-consistent spin–orbit total energies to within 5 µeV.

• MAE values computed before today are biased, by about 10% for FePt and by an amount that depends on the cell. Cached results are recomputed automatically.

How the route works

The route uses the magnetic force theorem through TB2J:

1. SCF: an ABACUS noncollinear SCF with spin–orbit coupling (SOC) switched off (soc_lambda = 0).

2. SOC step: a single ABACUS step with SOC switched on, on the same density.

3. Band energy: TB2J takes both Hamiltonians and evaluates the second-order SOC energy for each magnetization direction as a contour integral over the Green's function.

The MAE is the spread of those energies across directions.

Finding 1: the A100 did nothing

MAE ran on an A100 worker. The GPU sat at 0–4% utilization. Both ABACUS steps run in the localized (LCAO) basis with a CPU eigensolver, and they ran as one MPI rank with 8 threads. The other DFT routes had already moved away from that layout, because it is 5–6× slower on these machines. The only part that touched the GPU was TB2J's "GPU-optimized" band-energy step. Its GPU utilization stayed under 4%, for the reason covered in Finding 2.

Finding 2: the band-energy step was a Python loop

On L1₀ FePt, TB2J's band-energy step took 161 s on the A100 and 277 s on CPU. That was 60–85% of the whole job. Two things caused it.

• A per-energy loop. The code looped over 100 contour energies × 3 directions × 384 k-points. At every energy point it rebuilt the Green's function and re-rotated the SOC matrix, and it moved small matrices to the GPU one at a time. But G(z) = V·g(z)·V† comes from eigenvectors that are already known. Moving the SOC matrix into that eigenbasis once per k-point and direction reduces the trace at every energy to elementwise arithmetic: Tr(G·dH·G·dH) = Σ_ab g_a g_b M_ab M_ba, with M = V†·dH·V. The new kernel takes 2.4 s and matches the old loop to 10⁻¹⁶ eV on real FePt output.

• A Hamiltonian rebuilt 768 times. HamiltonIO's split-SOC model recomputes the entire real-space Hamiltonian every time it is read, and TB2J read it once per k-point, twice per run. Building it once was pure caching, and the energies came out bit-identical. That removed about 65 s per FePt job.

Every FePt run below except the last used the old TB2J loop. The band-energy step (lightest segment) dominates all of them, on the A100 as much as on CPUs:

Result: MAE runs on 16 CPUs

MAE now runs on a 16-core CPU worker, one MPI rank per core, with k-point pools and the genelpa solver. There is no GPU in the worker. Every configuration in a case had to reproduce the same MAE; they agreed to within about 10 µeV/f.u.

With the new TB2J kernel, the ABACUS SCF is most of the job, and it scales with MPI ranks and k-point pools:

Using the GPU properly, with the cusolver eigensolver, was the fastest option only for the 16-atom cell (284 s). It was 12% faster than 16 CPUs there, at 3× the cost:

Finding 3: the MAE depended on the cell

As a consistency check, we ran FePt as its 2-atom cell and as a 2×2×2 supercell with an exactly equivalent k-mesh (8×8×6 folds onto 4×4×3). Those must give the same MAE per formula unit. They didn't: 3.16 against 2.90 meV/f.u. Run-to-run scatter is about 10 µeV, so the gap was systematic.

The cause was a band cutoff in TB2J. It keeps band indices only up to the last band that dips below E_F + 5.1 eV at some k-point. The second-order SOC energy is a sum over virtual transitions into the empty bands, so the cutoff drops part of it. Because the cutoff counts bands, in a supercell (where eight primitive k-points fold onto each k) it lands at a different energy: 18 eV above E_F in the 2-atom cell, 8 eV in the supercell.

With the cutoff removed, both cells agree exactly. To find out which number is actually right, we added a reference that uses no perturbation theory: self-consistent ABACUS SCFs with full SOC, magnetized along 001 and along 100, with the MAE taken from the total-energy difference.

The contour isn't a factor: doubling the number of contour points moves the result by about 10 µeV. We also tried diagonalizing with SOC directly for each direction, which gives 3.55; that is not a valid reference. Each method's error against the self-consistent reference (the dataset table has every variant, including the contour checks):

The 5.1 eV cutoff is still used for exchange and Curie temperatures. Those depend on states near the Fermi level, and those results are unchanged.

MAE now respects ferrimagnetic seeds and +U

The MAE route used to start every site of an element from the same positive moment. It ignored initial_magmoms, magCIF moments, and the default that sets Mn/Cr antiparallel to Fe/Co/Ni. So a ferrimagnet such as GdCo₅ was quietly computed from its ferromagnetic state on every axis. See the gadolinium basis post

• Seeds: MAE now seeds exactly like the collinear routes and rotates each site's signed moment onto the trial axis. Results report the seed under magnetic_sublattice.

• +U: hubbard_u is now accepted on MAE. U is ramped in during the SCF. The one-step SOC run applies the full U to the SCF's converged occupation matrix and refuses to run if that matrix is missing. Both steps were checked on GdCo₅. ABACUS runs DFT+U in the noncollinear mode, and the SOC step's log shows it starting from the SCF's fully polarized Gd 4f occupation matrix, not from zero.

Checked on GdCo₅

GdCo₅ is a ferrimagnet: Gd couples antiparallel to Co, and its anisotropy comes from the Co sublattice. We relaxed it with U = 6.7 eV on the Gd 4f shell, starting from the ferrimagnetic seed, then ran MAE on the relaxed cell at the same settings (150 Ry, 8×8×8 k-mesh, all bands).

The easy axis is right, and K is about 40% low. A spin–orbit calculation of the moments at the same settings gives Co orbital moments of 0.12–0.14 μB, against about 0.25 μB measured in RCo₅ compounds. GGA with spin–orbit coupling has no orbital-polarization correction, so it underestimates Co orbital moments, and Co anisotropy tracks them. Treat MAE on Co-based magnets as a lower bound.

Before the band-window fix, the same cell gave 2.44 MJ/m³. The cutoff bias was +3% here, against +10% for FePt. The two basal-plane directions, which are equivalent by symmetry, differ by 0.05 meV/f.u. (4% of the MAE) at this k-mesh.

If you use MAE numbers from Ouro

• Recompute anything from before 2026-09-29. The band-cutoff bias has no single sign or size; it depends on how bands fold in your cell. For FePt it was +10%. The 3.45 meV/f.u. FePt figure in the Curie temperature and MAE cutoff post

Notes for anyone running TB2J MAE themselves

• The TB2J changes are in the vendored copy in the ouro-dft repo, under TB2J/MAE.py and TB2J/green.py. They add a kernel="eigen" option, the cached real-space Hamiltonian, and a TBGreen(emax=...) band-window parameter. The unit tests check that the MAE is the same in a supercell with the equivalent k-mesh; that check catches this class of problem immediately.

• If a dense contour and a converged k-mesh still leave your MAE different between a cell and its supercell, look at the band window before anything else.

• bench_mae.py in the same repo reruns every measurement in this post: hardware and process layouts, kernel against the old loop on the same ABACUS output, band-window variants, and the self-consistent SOC reference.