  Vibronic coupling has a dramatic influence over a large number of molecular
processes, ranging from photo-chemistry, to spin relaxation and electronic
transport. The simulation of vibronic coupling with multi-reference
wavefunction methods has been largely applied to organic compounds, and only
early efforts are available for open-shell systems such as transition metal and
lanthanide complexes. In this work, we derive a numerical strategy to
differentiate the molecular electronic Hamiltonian in the context of
multi-reference ab initio methods and inclusive of spin-orbit coupling effects.
We then provide a formulation of open quantum system dynamics able to predict
the time evolution of the electrons' density matrix under the influence of a
Markovian phonon bath up to fourth-order perturbation theory. We apply our
method to Co(II) and Dy(III) molecular complexes exhibiting long spin
relaxation times and successfully validate our strategy against the use of an
effective spin Hamiltonian. Our study shed light on the nature of vibronic
coupling, the importance of electronic excited states in spin relaxation, and
the need for high-level computational chemistry to quantify it.
