Standard implementations of non-relativistic excited-state calculations compute only one component of spin multiplets (i.e., Ms =0 triplets), however, matrix elements for all components are necessary for calculations of experimentally relevant spin-dependent quantities. To circumvent explicit calculations of all multiplet components, we employ Wigner–Eckart’s theorem. Applied to a reduced one-particle transition density matrix computed for a single multiplet component, Wigner–Eckart’s theorem generates all other spin–orbit matrix elements. In addition to computational efficiency, this approach also resolves the phase issue arising within Born–Oppenheimer’s separation of nuclear and electronic degrees of freedom. A general formalism and its application to the calculations of spin–orbit couplings using equation-of-motion coupled-cluster wave functions is presented. The two-electron contributions are included via the mean-field spin–orbit treatment. Intrinsic issues of constructing spin–orbit mean-field operators for open-shell references are discussed and a resolution is proposed. The method is benchmarked by using several radicals and diradicals. The merits of the approach are illustrated by a calculation of the barrier for spin inversion in a high-spin tris(pyrrolylmethyl)amine Fe(II) complex.