expressions for these quantities was recently formulated in the adiabtic-Fock basis for on-the-fly quantum dynamics simulations. However, these quantities are rarely available in standard electronic structure packages, including the derivative couplings d αβ ), due to the complexity of obtaining the analytical expression for excited state electronic structure method. The analytic derivative couplings have only recently been developed for NAMD simulations for common excited state methods like TD-DFT 238,239 over the last decade or so and implemented in only a few electronic structure or NAMD packages. Additionally, a recent work indicated that the explicit formulation of the derivative couplings may not be needed and can in fact be approximated very accurately only using the diagonal gradients and potential energies. Recently, Zhang and Tretiak implemented analytic nuclear gradients on the dipole and simulated the photoexcited dynamics of the stilbene molecule. In this work, the authors modified the NEXMD software package to include the pQED Hamiltonian (see Sec. 3) at the Jaynes-Cummings level with all proper gradients required for this Hamiltonian (i.e., without DSE and making the rotating wave approximation). Additionally, the gradients on the potential energy surfaces, non-adiabatic couplings, and dipole gradients were achieved analytically at the TD-AM1 255 level of theory in the collective electronic oscillator (CEO) framework. Most importantly, the nuclear gradient on the bare transition dipole between the ground and excited electronic states was computed as µ 0α = Tr[ μ X0α ] in the atomic orbital {o, v} basis and can be understood as,