  In conventional \textit{ab initio} methodologies, phonons are calculated by
solving equations of motion involving static interatomic force constants and
atomic masses. The Born-Oppenheimer approximation, where all electronic degrees
of freedom are assumed to adiabatically follow the nuclear dynamics, is also
adopted. This approach does not fully account for the effects of broken
time-reversal symmetry in systems with magnetic order. Recent attempts to
rectify this involve the inclusion of the velocity dependence of the
interatomic forces in the equations of motion, which accounts for time-reversal
symmetry breaking, and can result in chiral phonon modes with non-zero angular
momentum even at the zone center. However, since the energy ranges of phonons
and magnons typically overlap, the spins cannot be treated as adiabatically
following the lattice degrees of freedom. Instead, phonon and spins must be
treated on a similar footing. Focusing on zone-center modes, we propose a
method involving Hessian matrices and Berry curvature tensors in terms of both
phonon and spin degrees of freedom, and describe a first-principles methodology
for calculating these. We then solve Lagrange's equations of motion to
determine the energies and characters of the mixed excitations, allowing us to
quantify, for example, the energy splittings between chiral pairs of phonons in
some cases, and the degree of magnetically induced mixing between infrared and
Raman modes in others. The approach is general, and can be applied to determine
the adiabatic dynamics of any mixed set of slow variables.
