A diabatization protocol that combines the recently developed model space diabatization scheme and a mean-field treatment for spin-orbit coupling is presented. The protocol is highly efficient and easy to use. Its robustness is demonstrated in the calculations of the diabatic spin-orbit matrix elements of PH3 +, AsH3 +, and SbH3 +. It captures subtle geometry-dependence of the spin-orbit couplings and provides wave function information to explain the origin of the dependence.