Here, d Sj (R(t ) ) is the non-adiabatic coupling between polaritonic states, d R(t ) d t = R(t ) are the classical nuclear velocities, and ρ Sj (t ) = c * S (t )c j (t ) is the reduced density matrix in the polaritonic basis. Between nuclear time-steps, the electronic DOFs are updated via direct propagation of the time-dependent Schrodinger equation as,