where TR is the nuclear kinetic energy, U 0 ( R) represents the state-independent part of the potential, and Ve ( R) is the state-dependent part of the potential. Further, R represents a nuclear DOF, {|n } represents a set of diabatic electronic states, and V nm ( R) = n| Ve ( R)|m is the matrix element of Ve ( R) in this diabatic representation. The identity operator Î = N n=1 |n n| represents the identity in the electronic Hilbert space. In order to avoid an exponential numerical scaling with the number of DOFs of the system, different approximate methods have been developed. Trajectory-based quantumclassical methods are among the most successful ones as a) Electronic mail: dbossion@ur.rochester.edu b) Electronic mail: pengfei.huo@rochester.edu they scale linearly with the number of DOFs and allow for a simple numerical propagation scheme. One of the most popular trajectory-based approach is the surfacehopping method, where an ensemble of classical trajectories hop among electronic states upon non-adiabatic transitions, mimicking the wavepacket branching dynamics. The other widely used approach is the Ehrenfest trajectory dynamics where the nuclear DOFs feel a timedependent mean field potential from the quantum subsystem's time evolution.