In the current literature, molecular reaction dynamics calculations are applied probably using one of the following theories.
In all these methods, constructing a reasonable Hamiltonian is the first challenge. For example, a Hamiltonian for a Morse oscillator cannot be changed into algebraic form simply like a harmonic oscillator can with ladder operators. It requires a series of transformations because the position operator is exponentiated.This gets more difficult when you include an interaction term. A prefferred approach is to combine classical translation trajectories with algebraic quantum reaction dynamics for the molecular internal degrees of freedom. This is exactly where my research focuses, a simple analytic formulation of reaction dynamics. The current goals are explicitly:
To find a reaction coordinate the common way is to use a NewtonRaphson minimization procedure to find extrema in the Hessian matrix eigenvalues. Then, form a functional that takes steps away in the direction of steepest descent. As you'll see most of what I do is an improvised version of these methods, but I hope to formulate a general algorithm eventually. I will start by choosing a surface in which the saddle point is easily calculated.
With the parameters A=2 and B=C=5, we find a saddle point at {0.79324,0.79324} by taking a slice of the surface along x=y. From this slice, we obtain a skewed Morse function of which the minimum is found and then scaled back accordingly.
Once we find the saddle point we can use a steepest descent procedure to find the "intrinsic path", this is reaction coordinate. The saddle point here is the red circle. As you can see the surface symmetry gives the coordinate the same curvature on both the reactant and product sides. This red line becomes our translation coordinate in the construction of the complete "natural" coordinate system. The problem is that the red line is a numerical entity, a list of vectors. It can easily be cast to analytic form with spline interpolation of a few select points along this numerical trajectory. This results in parametric equation vectors of the form u(x,y)={x(t),y(t)} and v(x,y)={x(t),y(t)}. It is these formulae that will be used to calculate the dynamics given initial conditions of the interacting molecular species. 
Here is a view of the explicit code used to formulate this image: Mathematica Notebook
Here are two coupled equations, one classical, one quantum. Solving these equations reveals much of the collision detail. The difficult part is parametrizing the coordinates while simultaneously solving these equations of motion.
Here is an example of a trajectory that comes from initial relative position and velocity of our colliding species. The yellow line is the "classical" trajectory for a possible reactive collision as it has enough energy to make it over the saddle point. The little red circle now represents the origin of our natural coordinate system. The orange arrow points in the direction of the u coordinate and the green in the v. As we move steadily along u the yellow line deviates an amount v from the local minimum. With the quantum equation of motion we do not need this "classical" deviation from the reaction coordinate(red path). With the time dependent Lie parameters in the WeiNorman ansatz the local well shape becomes a function of time as you move along the "classical" translation. 
Here are two important slices taken from the above potential energy surface. The blue is along x=10 or y=10 since the surface is symmetric. The purple plot is the slice taken along x=y. The changing shape means changing eigenvalues.The algebraic method is seen here to be most advantageous. The changing well shape cooresponds to changing su(2) parameters. With the chain su(2)>so(2) we can see the corresponding quantum numbers,{N,m} are associated with the depth and width respectively. Note: Letting N go to infinity brings the simple harmonic oscillator eigenvalues back. Also, the V is for potential and the green arrow is v. 

This algebraic model for the spectroscopy of a given slice is uber powerful in statics but escalates in difficulty at the saddle point for dynamical calculations. Although, at this point it may be wise to consider the system to be coupled Morse oscillators. It must be noted that the purple potential well shape represents a counterintuitve potential which gives rise to energy eigenvalues of 3 bodies in anharmonic oscillation. The Frenet frame can be seen here as the blue vector, normal to the surface, completes the set of moving orthogonal coordinates.
Kinetic theory can be applied to a molecular gas where each molecule behaves like a hard sphere. Collisions involve simple cross sections and are used easily to find reaction rates.Here is the procedure
Taken from Hyperphysics