First-principles calculations were carried out using Density Functional Theory (DFT). All DFT calculations were performed within the VASP code , using the PBEsol generalised gradient approximation exchange-correlation functional . Projector augmented-wave pseudopotentials were used to model the ion cores , while treating the N 2s and 2p, and the Li 1s and 2s electrons as valence electrons. We used 4 × 4 × 4 Monkhorst-Pack k-point meshes to sample the Brillouin zone . k-point mesh sizes were reduced for relevant supercell calculations to preserve net k-point density. For geometry optimisations and force calculations, the electronic wave functions were expanded in a plane-wave basis with a kinetic-energy cutoff of 950 eV. The kinetic-energy cutoff was reduced to 400 eV for AIMD and the MLFF training. Geometry optimisations were performed using a quasi-Newton algorithm with a convergence setting of 0.005 eV Å-1 for the atomic forces. The resulting optimised structure of Li 3 N is shown in Figure . The calculated lattice parameters for this DFT-optimised structure compare well with previous computational and experimental data. More detailed information can be found in reference .