Coupling molecules to the quantized radiation field inside an optical cavity creates a set of new photonmatter hybrid states, which are commonly referred to as polaritons. These polariton states have hybridized curvatures from both the ground and the excited electronic states, which have been shown to facilitate new chemical reactivities. Thus, polariton chemistry provides a new paradigm for chemical transformations. Theoretical investigations play a crucial role in understanding the fundamental limit and basic principles in this emerging field, as these polariton chemical reactions often involve a rich dynamical interplay among the electronic, nuclear, and photonic degrees of freedom (DOFs). Accurately simulating polaritonic quantum dynamics remains a challenging task and is beyond the knowledge of photochemistry which does not include quantized photons, or quantum optics which does not have a well-defined theory to include the influence of nuclear vibrations. The mixed quantum-classical (MQC) approaches play an important role in studying the dynamics of molecules and simulating the non-adiabatic dynamics of the coupled electronic-nuclear DOFs. Two of the most commonly used MQC methods are Ehrenfest and fewest switches surface hopping (FSSH) approaches. Both approaches describe the electronic subsystem quantum mechanically, and treat the nuclear DOF classically. It is thus a natural idea for the theoretical chemistry community to extend these two approaches to investigate polariton chemistry by treating the electronic-photonic a) Electronic mail: pengfei.huo@rochester.edu DOF (or so-called polariton subsystem) quantum mechanically, and describing the nuclear DOF classically. Incorporating the description of the photon field into the MQC methods has become a basic problem of crucial importance to study polariton chemistry. The key ingredient in the MQC simulations of polariton dynamics is the expression of the nuclear gradient, which is a necessary ingredient for propagating the motion of nuclei. There has been large progress during the past few years along this direction. Kowalewski and Mukamel derived the expression of the derivative couplings using the Jaynes-Cummings (JC) model Hamiltonian where the rotating wave approximation (RWA) is assumed for the molecule-cavity coupling term. Groenhof and coworkers developed a multi-scale simulation approach combining Tavis-Cummings (TC) model using FSSH approach or Ehrenfest dynamics to simulate an ensemble of molecules coupled to a cavity. Zhang et al. extended the JC and TC models to include multiple molecular excited states, derived the corresponding nuclear gradient and performed MQC simulations for polariton dynamics. Fregoni et al. developed the MQC simulations with the JC-type model (that excludes the permanent dipole moment) to perform FSSH simulations of azobenzene photoisomerization in cavity. However, these previous simulations either used RWA, excluding the dipole-self energy (DSE) terms, or neglected permanent dipole moment, all of which could change the fundamental polariton chemistry and physics in the strong and ultra-strong coupling regime. Nuclear gradient expressions beyond the JC models has been derived under the full Configuration-Interaction (CI) expansion framework. Rigorous nuclear gradient expressions under the quantum electrodynamics Linear-Response time-dependent Density Functional Theory (QED-LR-TDDFT) framework using Pauli-Fierz (PF) type QED Hamiltonian. However, these gradient expressions lack a clear physical picture of light-matter interactions, as well as their clear connection to the more intuitive (but less accurate) gradient of the JC-type Hamiltonian.