MAE without the second SCF: the spin–orbit operator is geometry

Our magnetic anisotropy route computes K₁ from the second-order spin–orbit band energy: two ABACUS Hamiltonians, H₀ without spin–orbit coupling and H₁ with it, and a sum over the Brillouin zone of |⟨a|H₁ − H₀|b⟩|² (f_a − f_b)/(ε_a − ε_b). Until today H₀ came from a second self-consistent run, nspin=4 at soc_lambda 0, on top of the collinear SCF the route already does for its geometry gate. That second SCF was most of the route's DFT time: 106 s on L1₀ FePt against 15 s for everything else, 711 s on a 108-basis Fe–W cell.

It isn't needed. Both Hamiltonians are functions of the density alone, and the nspin=4 SCF converges to the same density as the collinear one (total energies agree to 10⁻¹⁰ Ry). So the route now takes two single nspin=4 steps on the collinear density, at soc_lambda 0 and 1, and sums as before. Same physics, same cache, same response; results carry h0_source: collinear_scf.

What the spin–orbit operator looks like

The detour that led here was a look at H₁ − H₀ itself across 141 cached MAE runs on 3d/4d/5d hosts. In the LCAO basis it is exactly what the pseudopotential says it should be:

• On-site: ξ L·S per atom and l, with ξ a per-element constant. Fe 3d 56.1 ± 0.15 meV across 134 jobs, Co 71.5, Pt 5d 518, W 285, Ta 225. The fit leaves 10⁻⁶–10⁻³ of the on-site norm, nothing in the s block, nothing between different l. No environment dependence at all: it is a table.

• Off-site: two-center terms from neighbor orbital tails at a heavy nucleus. Only 1% of ‖H₁ − H₀‖² for B/N/Al hosts and 5–7% for W/Ta, but the band sum is quadratic, so their cross term with the on-site part is first order in amplitude. Dropping them overshoots K₁ by 25% on FePt and by 1.5–3× on Fe–W and Fe–Ta, while 3d/sp hosts stay within 3%.

So the on-site table is a good approximation for 3d chemistry and a bad one for the 5d hosts we care most about, and the two-center part has to be computed. Since the whole operator is density-independent, that costs one single-step run, not an SCF.

Validation

Route (nspin=4 SCF) against the new recipe on cached structures, sum on the route's own mesh and smearing. K is the out-of-plane anisotropy in MJ/m³.

17 structures: max |Δ| 0.026 MJ/m³, median |Δ| 0.004. Route times are the cached nspin=4 SCFs on the 16-core MAE worker; the steps ran on 8 cores.

Each step is one diagonalization pass, 5–30 s on 8 cores. The collinear SCF they read is the one the geometry gate already made, so it costs the route nothing extra.

Limits

• DFT+U keeps the nspin=4 SCF: the collinear run's occupation matrix cannot seed the spinor one (h0_source: nspin4_scf).

• Rare earths are a limit of the theory, not the recipe. Second-order perturbation in spin–orbit coupling does not describe 4f shells.

• Lifting the nspin=2 Hamiltonian straight into spinor form, without the soc_lambda 0 step, is cheaper still but not good enough: the scalar-relativistic parts of the two formalisms differ at the 10⁻³ level, which moves small-K cells by up to 0.1 MJ/m³.

Why this matters beyond speed

H₀ is now a plain collinear Hamiltonian, which is what Hamiltonian-learning models (DeepH-style) are built to predict, and every cached collinear SCF that wrote H(R) is an anisotropy label for the price of two short steps. That is the path to a learned MAE screen that we could not get from public high-throughput datasets, where held-out-chemistry correlation never rose above 0.2.