As pointed out in previous sections (e.g., Sec. 4), performing trajectory-based ab initio polaritonic dynamics requires multiple properties of the molecular system, including the energy matrix, ( Ĥ ) αβ , transition dipole matrix, ( μ ) αβ , and the nuclear gradients of both ( R Ĥ ) αβ and ( R μ ) αβ . The gradient on the energy leads to the usual state gradients R E α as well as the first-order non-adiabatic couplings R d αβ (Eq. 4). These non-adiabatic couplings may be analytically available in some electronic structure software, such as QCHEM, 169 MOLPRO, 170;171 NEXMD, 155 etc. but not all, and largely depend on the level of electronic structure itself. Further, the nuclear gradients of the transition dipole matrix elements, to our knowledge, is not implemented in an analytical form in any common electronic structure software.