Two barriers, one tunnel: the hydrogen peroxide torsion, recomputed

Rebuilt July 2026. The original version of this post (“From Mock to Reality…”) documented moving a tunneling workflow from mock data to psi4, and its two lessons about modern psi4 usage survive below. Its physics did not. It called its 6.70 kcal/mol barrier the trans barrier and validated it against “the experimental trans barrier of ~7.4 kcal/mol” — but hydrogen peroxide’s big barrier is cis; the trans barrier is a factor of seven smaller. Its WKB and SCT columns were identical because the code aliased them, and its Eckart transmission was a step function — an implementation bug reported as a method limitation. This rewrite keeps the subject and the URL, recomputes everything with scripts small enough to show, and grades the results against the far-infrared spectroscopy of Hunt and co-workers¹ rather than against a misremembered number.

Hydrogen peroxide is the smallest molecule with a torsion, and its torsion is famously strange. The equilibrium structure is skew: the two O–H bonds sit at a dihedral angle near 112°, not planar. Planarity can fail in two ways — the cis structure (dihedral 0°, hydrogens eclipsed) and the trans structure (180°, hydrogens opposed) — and these are both saddle points, with famously unequal heights. Getting from one equivalent skew well to its mirror image is a chemical motion with no bond broken, a barrier small enough to tunnel through, and sixty years of precise spectroscopy to keep a calculation honest¹.

That last part is the working rule carried over from the fundamentals rewrite: every computed number here meets either a closed form, an independent numerical method, or an experimental value before it gets used. Three of those checks fail loudly in what follows, and each failure is the interesting part.

1. The potential, computed properly this time

The scan is psi4 1.11², MP2/cc-pVTZ with frozen core: a free optimization of the skew minimum, then constrained optimizations marching the H–O–O–H dihedral from 0° to 180° in 10° steps, each point seeded with the previous point’s geometry. Symmetry does the other half: \(V(-\varphi) = V(\varphi)\). Two practical lessons from the original post survive intact. Modern optking replaced fixed_dihedral with ranged_dihedral, and setting a tight range recovers fixed-value scans; and high-symmetry points (0°, 180°) will crash mid-scan with symmetry-projection errors unless the geometry forces symmetry c1 with no_reorient. The whole loop:

for phi in grid:                      # 0..180 in 10-degree steps, chained
    psi4.set_options({**BASE,        # basis, freeze_core, convergence
        "ranged_dihedral": f"3 1 2 4 {phi-0.01} {phi+0.01}"})
    geom_from_zmat(phi, *seed)        # seeded with the previous geometry
    E = psi4.optimize("mp2")
    seed = internals_of_current()     # r_OO, r_OH, angle for the next point

The optimized minimum lands at \(\varphi_{\mathrm{eq}} = 114.3°\), with \(r_{\mathrm{OO}} = 1.450\) Å, \(r_{\mathrm{OH}} = 0.964\) Å — a couple of degrees and a few thousandths of an ångström from the CCSD(T) reference values³. An eighth-order cosine series fits all nineteen scan energies with a maximum residual of \(0.07\ \mathrm{cm^{-1}}\), and MP2 frequency calculations at both planar structures find exactly one imaginary mode each — \(307.9i\ \mathrm{cm^{-1}}\) at trans, \(611.1i\ \mathrm{cm^{-1}}\) at cis — confirming both are true first-order saddle points of the torsion.

Figure 1. The computed torsional potential of H\(_2\)O\(_2\) (black; the cis walls at 0° and 360° continue up to \(2670\ \mathrm{cm^{-1}}\), off the top of the frame) with the four lowest torsional levels from §2 drawn between their classical turning points — symmetric states in blue, antisymmetric in red. The ground pair hugs the wells and its two lines nearly coincide (\(13\ \mathrm{cm^{-1}}\) apart); the first excited pair straddles the trans barrier top (dashed) and is split by \(117\ \mathrm{cm^{-1}}\).

Table 1. Stationary points of the torsion against the empirical hindering potential fitted to the far-infrared spectrum. Koput’s CCSD(T)/cc-pVQZ surface puts the trans barrier near \(375\ \mathrm{cm^{-1}}\)³ — our MP2 value is serviceable, a little low.

This table is where the original post went off the rails, so it is worth being explicit. The big barrier is cis — eclipsing the hydrogens costs lone-pair repulsion — and the small barrier is trans. The old post reported one barrier of 6.70 kcal/mol, called it trans, and matched it against “~7.4 kcal/mol”, a number that belongs to the cis barrier. The agreement it celebrated was a coincidence of two errors. Worse, the mislabeling propagated: every tunneling quantity downstream was computed as if the torsion had to cross a 7 kcal/mol wall, when the road from one well to the other over trans costs one-seventh of that.

2. The observable: a splitting, not a rate

What does tunneling actually do to this molecule? H\(_2\)O\(_2\) in its ground state does not “react” — it sits in one skew well, and the well’s mirror image sits \(131°\) of dihedral away. Tunneling connects them, and the signature is not a rate constant but a level splitting: each torsional state of an isolated well becomes a symmetric/antisymmetric pair, separated by an energy that measures how often the molecule leaks through the wall. Hunt’s far-infrared spectrum resolved exactly this ladder in 1965¹.

The test is therefore parameter-free: solve the one-dimensional torsional Schrödinger equation on our computed potential,

\[ -\frac{\hbar^2}{2I}\frac{d^2\psi}{d\varphi^2} + V(\varphi)\,\psi = E\,\psi, \qquad \psi(\varphi + 2\pi) = \psi(\varphi), \]

and compare the resulting splittings with the measured ones. The effective moment of inertia \(I\) comes from the scan geometries themselves (two O–H rotors counter-rotating about the O–O axis; \(I = 0.456\) amu Ų at the minimum, varying by only ±2% along the path). In a periodic plane-wave basis the Hamiltonian is a small dense matrix — the same eigenvalue-problem machinery as every post in this series, eight lines of NumPy:

m = np.arange(-80, 81)                       # plane waves e^{im phi}
H = np.diag(B*m.astype(float)**2 + c[0])     # B = hbar^2/2I in cm^-1
for n in range(1, 9):                        # V as a cosine series
    H += 0.5*c[n]*(np.eye(161, k=n) + np.eye(161, k=-n))
levels = np.linalg.eigvalsh(H)

Table 2. Torsional energy levels relative to the \(0^+\) ground state, from the plane-wave solve on the MP2/cc-pVTZ potential, against the far-infrared values. Ground-state tunneling splitting: computed \(13.0\ \mathrm{cm^{-1}}\), observed \(11.4\ \mathrm{cm^{-1}}\). The first-excited splitting comes out \(117.4\ \mathrm{cm^{-1}}\) against an observed \(116.5\ \mathrm{cm^{-1}}\).

For a rigid one-dimensional model with no adjustable parameters, this is a startlingly good scorecard, and the errors have honest signs: our trans barrier is \(31\ \mathrm{cm^{-1}}\) too low, so the molecule tunnels a little too easily and the ground splitting comes out \(14\%\) large. Varying \(I\) over its full along-path range moves the splitting only between \(12.8\) and \(13.3\ \mathrm{cm^{-1}}\) — the barrier, not the moment of inertia, owns the error budget.

Figure 1 shows why the two splittings differ by an order of magnitude. The ground pair sits at \(160\ \mathrm{cm^{-1}}\), less than halfway up a \(355\ \mathrm{cm^{-1}}\) wall — genuine deep tunneling, and the pair’s two lines are nearly indistinguishable. The first excited pair straddles the barrier top, where “tunneling” is shading into “barely hindered rotation,” and the pair splits wide open. One more encounter with an old friend: the tall cis wall is doing real work in this spectrum — it is what keeps the levels from being free-rotor states — but as a tunneling pathway it is irrelevant, as §3 quantifies.

3. Transmission and the honest kappa

The original post’s headline was a tunneling correction \(\kappa = 3.64\) at 300 K, computed — recall — through the wrong barrier with two of its three methods secretly the same method. Redo it properly. The semiclassical WKB transmission through the computed trans barrier (Kemble form, barrier penetration integral evaluated on the fitted potential with the along-path \(I(\varphi)\)) and the analytic Eckart transmission⁴ matched to the same barrier height and the psi4 imaginary frequency:

Figure 2. Transmission through the trans barrier on a log scale: WKB on the fitted potential (blue) and the analytic Eckart barrier matched to the computed height and imaginary frequency (red). The dashed line marks the barrier top. Both curves are smooth through \(E = V_b\) with \(T(V_b) \approx 0.5\)–\(0.6\) — not a step function. Where the old post’s Eckart jumped from 0 to 1 at the barrier top, the analytic result spreads that jump over roughly \(\hbar\omega^\ddagger \approx 300\ \mathrm{cm^{-1}}\).

Two verification notes, in the spirit of this series. The Eckart curve here is not trusted because the formula looks right: the script also computes transmission through the same \(\mathrm{sech}^2\) barrier by direct numerical integration of the Schrödinger equation (Numerov, matched to plane waves on both sides), and the analytic and numerical values agree to five decimal places across the energy grid. That check is what caught a factor-of-two error in my own first implementation of the barrier width — the same class of bug as the old post’s step function, caught the way such bugs have to be caught, by an independent method rather than by eyeballing plausibility.

With \(T(E)\) in hand, the thermal tunneling correction is the Boltzmann average \(\kappa(T) = e^{V_b/k_BT}\!\int_0^\infty T(E)\, e^{-E/k_BT}\,dE\,/\,k_BT\)⁵:

Figure 3. The tunneling correction \(\kappa(T)\) for crossing the trans barrier: WKB on the computed potential, analytic Eckart, and the Wigner low-order estimate \(1 + (\hbar\omega^\ddagger/k_BT)^2/24\). Three genuinely different methods now agree within ~10% — compare the old post, where two “different” methods agreed to six digits because they were the same code.

Table 3. \(\kappa(T)\) for the trans path. At room temperature, tunneling is a 7–15% correction, not the old post’s factor of 3.64; even at 100 K it is a factor of ~2, not ~47. A 1D transition-state estimate with the computed torsional frequency as prefactor puts the classical interwell rate near \(2\times 10^{12}\ \mathrm{s^{-1}}\) at 300 K — the torsion flips on the picosecond scale, four orders of magnitude faster than the old post’s \(1.3\times 10^{8}\ \mathrm{s^{-1}}\), which was computed against the wrong barrier.

And the cis pathway? Its \(\kappa\) is formally enormous at low temperature (tunneling deep under a tall barrier beats the classical route over it), but multiplying by the Boltzmann factor of a \(2670\ \mathrm{cm^{-1}}\) wall settles the question: the tunneling-corrected cis flux is \(10^{-12}\) of the trans flux at 100 K and \(2\times10^{-5}\) of it at 300 K. Every microscopic reversal of this molecule’s handedness goes through trans. The tall barrier shapes the spectrum; the short one carries the traffic.

4. The autopsy, and the limits

What the original post got wrong compresses to three entries, each with a general lesson. The barrier was validated against the wrong experiment — cis and trans swapped — and a validation against the wrong number is worse than none, because it converts an error into a certificate. Two methods agreed because they were one method — the SCT column was the WKB column under a second name; agreement between methods is evidence only when the methods are independent. A broken implementation was reported as a method limitation — the step-function “Eckart” — when a ten-line numerical scattering check would have convicted the code instead of Carl Eckart.

The rebuilt numbers carry their own stated limits. This is a rigid 1D model: one torsional coordinate, a moment of inertia frozen per-geometry, no zero-point adjustment of the other five modes along the path, and an MP2 barrier that sits \(31\ \mathrm{cm^{-1}}\) below the empirical one. Those approximations show up honestly as the 14% overshoot in the ground splitting. The systematic fix — letting the path curve through all coordinates and treating deep tunneling variationally — is instanton territory^5,6, a genuinely different machine. The point of this rewrite is narrower: with the right barrier labeled by the right name, a laptop-scale scan plus an eigenvalue problem reproduces sixty-year-old far-infrared spectroscopy to a few percent, and every claim in the post is attached to the check that would catch it being wrong.

References