Recently, various theoretical methods have been developed or extended to directly simulate the polariton dynamics. These include the full quantum dynamics simulations, mixed-quantum-classical (MQC) dynamics, and non-adiabatic dynamics based on the mapping formalism. Among them, the MQC dynamics methods describe the electronic-photonic DOFs quantum mechanically and treat the nuclear DOFs classically, hence well balance the computational cost and accuracy of the dynamics. Thus, the MQC methods, including the Ehrenfest and trajectory surface hopping (TSH) methods, have been widely used in the nonadiabatic polariton dynamics recently. In the propagation of the polariton dynamics with the MQC methods, besides the energies of the electron-photon hybrid states, we need to derive the nuclear gradients and the couplings between these states, where the derivatives of molecular dipoles (including permanent dipoles and transition dipoles) are the key ingredients. For some model systems with well-defined diabatic electronic states, the dipoles of/between these diabatic states can be set to constants. For model systems with adiabatic electronic states, the dipoles of/between these adiabatic states can be calculated through the discrete variable representation (DVR). Evaluating the derivatives of molecular dipoles remains a theoretical bottleneck for simulating ab initio polariton quantum dynamics. These derivatives are neither readily available for most electronic structure methods nor computationally cheap to obtain.