How does one find the probability of tunneling along a reaction coordinate computationally? The idea for this question comes from Jan's comment over here. I realized I don't know how I would compute the probability of tunneling in the process of a reaction.

When I say "compute" I don't mean something rather trivial like finding the probability of a quantum harmonic oscillator being past the classical turning points which can just be done analytically.

For instance, if I were studying some reaction computationally (to estimate rate constants or something which depends on tunneling) and constructed a potential energy surface for a bond-breaking process, let's say the transfer of a hydrogen, how do I account for the possibility of tunneling through the barrier along this reaction coordinate?

I could speculate, but I imagine I would be wrong. How would one do this using common computational software and techniques?
If the reaction coordinate $\xi$ can be expressed as the change in the relative position between two groups of the molecule (e.g. for instance as a change in a bond length so that $\xi$ is given in units of length) the problem may be treated as one dimensional and one may use computational software to calculate the potential energy $V(\xi)$ along this path. You may then proceed solving the time-independent Schroedinger equation using numerical techniques such as Numerov's integration method. However, in many cases WKB theory provides an accurate estimate. 

In terms of WKB, the tunneling probability is given by
$$
T=e^{-2G}
$$
where $G$ is the Gamow factor given by
$$
G=\frac{2}{\hbar}\int_{\xi_a}^{\xi_b}{ \sqrt{2\mu_\text{red}\left [V(\xi) - E \right ]}}d\xi
$$
where $\mu_\text{red}$ is the reduced mass along the reaction coordinate and $E$ is the energy of the incoming particles. The integral runs over the clasically forbidden region from $\xi_a$ to $\xi_b$. The Gamow factor thus represents the square root of the integral over the "reaction barrier" multiplied with the effective mass that tunnels. Of course, you should take zero point energy into account and in the case of light reactants also the vibrational and rotational structure of the products as this may lead to resonances.