We use Tully's FSSH algorithm to perform surface hopping simulations for polariton dynamics. Note that a similar simulation has been performed recently, and our focus here is to test the importance of the new gradient expressions. The equation of motion for electronic wavefunction in FSSH is described in Eq. 19, which is same as Ehrenfest Dynamics. The main difference is how the nuclear force is treated. In FSSH, the nuclear force is proposed to be contributed from only one specific eigenstate of V (so-called the active state), i.e., Eq. 42. The active state index will be determined at every single nuclear propagation step. According to the "fewest switches" algorithm, the probability of switching (probability flux) from the active polariton state |I to any other polariton state |J during the time interval between t and t + δt is