resulting states ( 6 H15/2 multiplet) are projected onto a crystal field (CF) Hamiltonian acting in the (2J+1)-dimensional |mJ⟩ basis, using our angmom_suite code (version 1.12.0). The spinphonon coupling parameters for each phonon mode, defined by band index j and wavevector q, are the derivatives of the CF parameters along the phonon normal mode vectors. To compute these, we first obtain the derivatives of the CF parameters with respect to Cartesian atomic coordinates using an analytic linear vibronic coupling (LVC) model, and then convert to the normal mode basis using the linear combination of atomic displacements specified by the mode displacement vector. This is done for all phonon modes, at all q-points in the sampling mesh.