The change in electronic state from one spin multiplicity to another, known as intersystem crossing, occurs in molecules via the relativistic phenomenon of spin-orbit coupling. Current means of estimating intersystem crossing rates rely on the perturbative evaluation of spin-orbit coupling effects. This perturbative approach, valid in lighter atoms where spin-orbit coupling is weaker, is expected to break down for heavier elements where relativistic effects become dominant. Methods which incorporate spin-orbit effects variationally, such as the exact-two-component (X2C) method, will be necessary to treat this strong-coupling regime. We present a novel procedure which produces a diabatic basis of spin-pure electronic states coupled by spin-orbit terms, generated from fully variational relativistic calculations. This method is implemented within X2C using time-dependent density-functional theory and is compared to results from a perturbative relativistic study in the weak spin-orbit coupling regime. Additional calculations on a more strongly spin-orbit-coupled [UO2Cl4]2− complex further illustrate the strengths of this method. This procedure will be valuable in the estimation of intersystem crossing rates within strongly spin-coupled species.