Imagine preparing a molecule with a short pulse of light. Some molecules return to their original structure; others rearrange or break a bond. A barrier on the ground-state energy surface cannot, by itself, tell us which outcome wins. The molecule was prepared elsewhere, and its electrons and nuclei can change their motion together before a product is formed.

Part II asked how a landscape becomes a reaction rate. Here the question changes: what if one landscape is not enough? We follow a major line of exploration from trajectory surface hopping in the early 1970s to later mixed quantum–classical and wavepacket methods. This series began with hydrogen transfer, but its subject is now the broader history of calculating chemical reactions. There is no chosen method waiting at the end of the story.

1. A new question: where does the population go?

For a thermally activated reaction on one electronic state, a useful target is a rate constant. For an electronically excited molecule, equally natural targets are the lifetime of an excited population, the fractions reaching different products, and the energy carried by each outgoing fragment. Electron transfer and changes of spin introduce related questions, although spin-changing dynamics also require the appropriate spin–orbit interactions.

These questions predate the 1970s: the Born–Oppenheimer framework and the Landau–Zener crossing models already supplied important foundations. The historical step emphasized here is making coupled electronic and nuclear motion into a practical molecular simulation. In 1971, Tully and Preston applied trajectory surface hopping to the reaction of a proton with a deuterium molecule. In 1990, Tully introduced the fewest-switches algorithm used in the example below. These are distinct milestones, not two dates for the same algorithm. See the original 1971 paper and 1990 paper.

The central compromise is easy to state: retain quantum amplitudes for electronic states, while describing the much larger set of nuclear coordinates with trajectories. The difficult part is making the two descriptions influence each other consistently.

2. What changes when one surface is not enough?

MQC 0 develops the Born–Huang expansion and derives the coupling terms. The essential point for this historical question is that nuclear motion differentiates the geometry-dependent electronic basis as well as the nuclear wavefunction. Keeping several electronic states therefore produces coupled nuclear channels, rather than independent trajectories on unrelated energy curves.

For nondegenerate states, the derivative-coupling expression summarizes where the single-surface picture can become fragile:

\[ d_{ij}(x)= \frac{\langle\phi_i|\partial_x\hat H_e|\phi_j\rangle} {E_j-E_i},\qquad i\ne j. \]

A small gap can amplify coupling, but the numerator and the speed and direction of nuclear motion still matter. At an exact degeneracy this formula is not valid. Most importantly, a crossing geometry alone does not supply a product yield or a population lifetime. Those require a dynamical calculation of where population goes.

3. Two representations of the same problem

A small two-state model keeps the distinction concrete. In a fixed diabatic basis, write the potential as

\[ V(x)=\begin{pmatrix}z(x)&u(x)\\u(x)&-z(x)\end{pmatrix}. \]

The off-diagonal potential \(u\) couples the two basis states. Diagonalization gives the adiabatic curves. Indeed,

\[ \begin{aligned} \det[V-EI]&=E^2-z^2-u^2=0,\\ E_\pm(x)&=\pm\sqrt{z(x)^2+u(x)^2}. \end{aligned} \]

When the diabatic diagonal elements cross at \(z=0\), the adiabatic gap is \(2|u|\). The avoided crossing and the diabatic crossing are thus not contradictory pictures. One uses diagonal electronic energies and derivative couplings; the other uses a coupled potential matrix. They describe the same model when every transformed term is retained.

For completeness, choose a continuous mixing angle \(\theta=\operatorname{atan2}(u,z)\), with lower and upper eigenvectors \((-\sin\frac\theta2,\cos\frac\theta2)\) and \((\cos\frac\theta2,\sin\frac\theta2)\). Differentiation yields

\[ d_{-+}=\frac{\theta'}{2} =\frac{zu'-uz'}{2(z^2+u^2)}. \]

Reversing an eigenvector's sign also changes the corresponding coupling and amplitude signs; observables remain unchanged if the convention is consistent. Globally eliminating derivative couplings is not always possible for a multidimensional molecule. Our fixed two-state diabatic Hamiltonian is a deliberately simple model, not a claim that every molecular diabatization is exact.

4. One average force, or different outgoing branches?

The moving-basis electronic equation already tells us how amplitudes evolve along a prescribed nuclear trajectory. The unresolved approximation is the feedback: which force should those amplitudes exert on the nuclei?

Ehrenfest dynamics uses the expectation of the electronic force operator. This includes coherence contributions; it is not generally just a population-weighted average of surface gradients. The force is smooth, but one mean-field trajectory cannot split into distinct outgoing nuclear packets. If one branch reflects and another transmits, the average trajectory need not represent either outcome.

Surface hopping chooses a different closure: each trajectory follows one active surface, while an ensemble can separate into branches. Electronic amplitudes and active-surface labels are different bookkeeping layers. This is a practical way to represent branching, not an exact rewriting of the coupled wavefunction.

5. The hopping rule—and the choices behind the benchmark

The MQC I derivation obtains the fewest-switches rule from outgoing electronic population flow. In the real adiabatic convention used here, recall

\[ \begin{aligned} J_{a\to b}&=2\operatorname{Re}(c_a^*c_b\dot x d_{ab}),\\ g_{a\to b}&=\max\left(0, \frac{J_{a\to b}\Delta t}{|c_a|^2}\right). \end{aligned} \]

The new issue is what this rule predicts for an experimentally motivated outcome, not the algebra of deriving it again. Active-state fractions need not exactly track averaged electronic populations once trajectories diverge, and frustrated hops create an additional inconsistency.

To make the comparison reproducible, our calculation uses the standard active-surface force, no decoherence correction, and energy-conserving momentum rescaling after an accepted hop. An energetically impossible upward hop is rejected with momentum unchanged. Other frustrated-hop prescriptions can change results; the rule is part of the model being tested. Time steps must resolve the population flow rather than routinely clipping probabilities above one.

The computational loop is: sample the incoming ensemble; advance nuclei and electronic amplitudes; draw hops and rescale accepted momenta; continue until channels separate; count reflected and transmitted outcomes on each surface. More trajectories reduce sampling noise, but cannot remove the physical approximations.

6. A calculation small enough to audit

We use Tully's simple avoided-crossing model, in atomic units, with \(M=2000\):

\[ \begin{aligned} z(x)&=A\operatorname{sgn}(x)(1-e^{-B|x|}),\\ u(x)&=Ce^{-Dx^2},\\ A&=0.01,\quad B=1.6,\\ C&=0.005,\quad D=1. \end{aligned} \]

The nuclei arrive from the left on the lower electronic state. We compare three descriptions: classical motion confined to that lower surface, fewest-switches surface hopping without a decoherence correction, and a two-component quantum wavepacket. The potential, mass, and incoming distribution are the same in all three. These are new calculations of the standard model, not digitized data from the original paper.

Match the preparation before comparing the answers

The incident quantum packet has mean position \(x_0=-12\), width \(\sigma=1\), and mean momentum \(p_0\):

\[ \chi(x,0)=(2\pi\sigma^2)^{-1/4} e^{-(x-x_0)^2/(4\sigma^2)} e^{ip_0(x-x_0)/\hbar}. \]

Squaring the envelope gives a position variance \(\sigma^2\). Fourier transformation gives momentum variance \(\hbar^2/(4\sigma^2)\). Its positive Gaussian Wigner function therefore allows trajectory initialization with independent samples

\[ \begin{aligned} x&\sim\mathcal N(x_0,\sigma^2),\\ p&\sim\mathcal N\!\left(p_0,\frac{\hbar^2}{4\sigma^2}\right). \end{aligned} \]

The packet starts in diabatic component 1, which is indistinguishable from the lower adiabatic state this far from the coupling region. The benchmark uses \(p_0=6,12,20,30\), with 4,000 trajectories per momentum and method. A finite-width packet averages over incident energies. Comparing it with one strictly monochromatic trajectory would mix a preparation difference into a method comparison.

What is “exact” here?

For the quantum reference, propagate the two-component wavefunction under the same Hamiltonian:

\[ \hat H=-\frac{\hbar^2}{2M}\partial_x^2 I+V(x). \]

A Fourier grid, also called a Fourier DVR, represents position and momentum on dual grids. The kinetic step is diagonal in momentum; the potential step is a local two-by-two matrix in position. Symmetric splitting follows by matching the short-time exponential expansion through second order:

\[ \begin{aligned} e^{-i\hat H\Delta t/\hbar} &=e^{-iV\Delta t/(2\hbar)} e^{-i\hat T\Delta t/\hbar} \\&\quad\times e^{-iV\Delta t/(2\hbar)} +O(\Delta t^3). \end{aligned} \]

There is no hopping rule in this propagation. Nuclear reflection, splitting, and coherent evolution all follow from the coupled Hamiltonian. The splitting has a finite time-step error, and the grid has a finite spatial error, so “exact” means a numerically converged result for this one-coordinate, two-state model. It does not mean exact chemistry or error-free arithmetic.

The reported reference uses 4,096 points on \([-64,64)\) and \(\Delta t=0.25\). We compare it with 2,048 points and \(\Delta t=0.5\), and with a larger box and later observation. The script also records probability near the box edges, unseparated central probability, and norm conservation. No absorber is used; observations must precede significant periodic-boundary return.

Count the same outgoing channels

Let \(r_V=\sqrt{z^2+u^2}\). Since \(V^2=r_V^2I\), its adiabatic projectors are

\[ P_\pm(x)=\frac12\left[I\pm\frac{V(x)}{r_V(x)}\right]. \]

This follows by applying the expression to an eigenvector with eigenvalue \(\pm r_V\): the matching projector returns one and the other zero. Once the packets have separated, integrate on either side:

\[ \begin{aligned} T_\pm&=\int_{x>4}\Psi^\dagger P_\pm\Psi\,dx,\\ R_\pm&=\int_{x<-4}\Psi^\dagger P_\pm\Psi\,dx. \end{aligned} \]

For trajectories, the corresponding quantities are fractions exiting to the right or left on the specified active surface. The electronic coefficient \(|c_+|^2\) on one hopping trajectory is not itself an outgoing-channel count.

Two-state potentials and quantum versus surface-hopping scattering probabilities
Left: the two diabatic diagonal potentials (dashed gray) and adiabatic eigenvalues (solid navy and rust), in hartree versus position in bohr. Middle and right: upper-state transmission and lower-state reflection versus mean incident momentum, in atomic units. Navy circles are quantum-grid results; rust squares are FSSH; gray triangles are trajectories restricted to the lower surface. Trajectory bars are 95% Wilson binomial sampling intervals, not total error bounds. The reflection panel uses an expanded vertical scale; gray and rust results there overlap. Lines connect calculated points as guides to the eye.
Mean \(p_0\) (a.u.)Quantum \(T_+\)FSSH \(T_+\)Quantum \(R_-\)FSSH \(R_-\)
6\(<10^{-8}\)0 observed0.022350.00175
120.21866\(0.2243\pm0.0066\)\(8.1\times10^{-7}\)0 observed
200.49299\(0.5088\pm0.0079\)\(8.6\times10^{-9}\)0 observed
300.71538\(0.7195\pm0.0071\)\(2.6\times10^{-10}\)0 observed

The table's \(\pm\) values are one sampling standard error, \(\sqrt{\hat P(1-\hat P)/N}\), obtained from the variance of a binomial count divided by \(N^2\). They are narrower than the figure's 95% intervals. Zero observed events does not prove a zero underlying probability: with 4,000 trials, the Wilson upper limit is about \(9.6\times10^{-4}\). Tiny quantum reflection values should likewise be read as negligible at this benchmark's accuracy, not as relatively converged predictions to every printed digit.

At \(p_0=12\) and 30, FSSH captures upper-state transmission to within the sampled uncertainty. At \(p_0=20\), the difference is about 0.016, roughly two standard errors; halving the time step gives 0.5025, so this small ensemble does not cleanly resolve a systematic bias there. The low-momentum reflection is a clearer failure: at \(p_0=6\), the quantum result is 2.23%, whereas both trajectory methods give 0.175%. Adding electronic hops has not supplied full quantum nuclear dynamics.

Why does upper-state transmission grow with speed? On the far left, diabatic state 1 is the lower adiabatic state; on the far right, it is the upper one. A fast passage can retain much of its original diabatic character while changing its adiabatic label. “An electronic transition” therefore needs a specified basis; it does not always mean a sudden change of chemical charge localization.

The largest absolute channel-probability change in the grid/time-step check is \(6.7\times10^{-7}\); the larger-box/later-time check changes a channel by less than \(1.2\times10^{-7}\). At most \(2.0\times10^{-6}\) of the quantum probability remains in the central region. Norm drift is below \(10^{-11}\), all trajectory exit counts are resolved, and the largest trajectory energy drift is below \(4.1\times10^{-7}\) hartree. These checks support the displayed probability precision, not a universal guarantee about either method.

7. What this comparison can—and cannot—establish

The one-surface model has no mechanism for populating the upper surface: its \(T_+\) is identically zero. Allowing electronic transitions opens a channel that was missing structurally, not merely underestimated through a slightly inaccurate barrier. Agreement with the coupled quantum result tests whether the hopping dynamics allocates population sensibly for this preparation.

It does not establish that surface hopping is universally accurate. Independent classical trajectories do not retain the relative nuclear phases needed for interference between separated packets. Their electronic amplitudes can also stay coherent after the associated nuclear branches have separated, motivating decoherence treatments. Frustrated-hop choices, classically forbidden nuclear motion, repeated passages, and electronic-structure errors can all matter in other regimes. Increasing the trajectory count only addresses statistical uncertainty.

Most importantly, the plotted numbers are probabilities, not rate constants. “Upper” is an electronic-state label and “transmitted” is a model scattering label; neither automatically denotes a chemically distinct product. To calculate an experimental reaction yield, a molecular application must define product structures and prepare the same initial ensemble as the experiment.

For an isolated excited population following a single exponential law, the bridge to a lifetime is simple. If \(dP_e/dt=-kP_e\), separating variables and integrating gives

\[ \begin{aligned} \ln\frac{P_e(t)}{P_e(0)}&=-kt,\\ P_e(t)&=P_e(0)e^{-kt},\qquad \tau=1/k. \end{aligned} \]

But coherent transfer, back transfer, or multiple channels need not produce an exponential. Nor is a measured transient absorption signal generally equal to a state population: transition strengths, nuclear geometry, pulse characteristics, and the instrument response enter the forward model. The observable-first lesson of Part I still applies.

8. An important exploration, not a final answer

Surface hopping is one answer to a cost problem, not the only answer to multistate dynamics. Another line keeps an explicit nuclear wavefunction while seeking a more economical representation. Meyer, Manthe, and Cederbaum's 1990 MCTDH work used time-dependent configurations to address multidimensional propagation. Multiple-spawning methods developed a different strategy: moving nuclear basis functions with additional branches generated when needed; Ben-Nun, Quenneville, and Martínez's 2000 article connects that approach to first-principles photochemistry. These developments pursued different balances between quantum detail and molecular scale.

Part IV turns to free energy and rare-event sampling: even with a usable Hamiltonian, how do we sample configurations that ordinary trajectories almost never see? After that, we will explore work on tunneling and quantum nuclear reaction dynamics. These strands overlap historically—Part IV returns to developments of the 1970s—and neither is simply a replacement for nonadiabatic dynamics.

Reproduce and inspect

python assets/code/reaction-dynamics/nonadiabatic_scattering.py

The script defines the Hamiltonian, initial ensemble, random seed, propagators, exit rules, error checks, and figures. More specialized implementation notes are available in the separate FSSH and Ehrenfest articles; their earlier benchmark setups are not interchangeable with the matched-ensemble calculation here.

Historical sources

Chemical Reaction Computation series map

  1. Part I: What Can We Actually Observe? — the series begins with hydrogen-transfer experiments.
  2. Part II: From a Barrier to a Rate — Arrhenius, TST, and a quantum rate benchmark.
  3. Part III: When One Electronic State Is Not Enough — multistate reactions and nonadiabatic dynamics.
  4. Part IV: Free Energy and the Events We Rarely See — equilibrium coverage, umbrella sampling, and the missing dynamical clock.
  5. Part V: From Barrier Penetration to Imaginary-Time Instantons — a thermal saddle, its ring representation, and an exact scattering benchmark.
  6. Part VI: Instanton Rates Across Crossover — Richardson, Lawrence, and the step from a saddle to a rate.
  7. Part VII: Multidimensional Instantons — a coupled 2D path solve, fluctuations, and numerical rate checks.
  8. Part VIII: Instantons and RPMD — statistical weight, recrossing, and a quantum-checked rate benchmark.