Count trajectories and measure the error

Benchmark Nonadiabatic Reaction Methods

Can more trajectories remove an approximation’s error in a reaction probability? Compare all three methods with the same Hamiltonian, initial packet and exit windows. Finer grids test the quantum reference, and more trajectories reduce sampling noise. A stable difference after these checks measures approximation error for that observable.

Compare channel probabilities, sampling error and computation time for quantum, Ehrenfest and fewest-switches dynamics. Save each calculation, change the ensemble size and see which differences survive.

Model assumptions

The reference is a converged numerical solution of a fixed two-state Tully model, not experimental chemistry data. Sampling errors describe independent Wigner preparations and hopping draws. Time-step and grid errors are separate. Propagation costs are measured on your device and exclude snapshots and rendering; they are not a universal ranking of the algorithms.

t=0

One nuclear coordinate. Height shows electronic energy; ribbon thickness shows density on a separate scale. Cyan is lower, coral is upper.

At the common final time

Separate sampling scatter from method bias

The exit windows begin four bohr from the crossing. Channel estimates include probability absorbed far beyond its exit; surviving and absorbed parts are saved separately. Probability between the windows is unresolved.

Quantum upper exit

0

Hopping minus quantum

0

Hopping standard error

0

Quantum unresolved

0

Four exit probabilities

Bars share the same probability scale. Thin whiskers show one sampling standard error for the trajectory estimates.

Saved estimates and their sampling error

Save runs with the same physical preparation to compare sample sizes. The line is the quantum value; each gold point is a saved hopping estimate.

Norm, energy, flux and measured cost

Probability still on the grid

0

Probability absorbed outside the exits

0

Norm plus absorbed

0

Quantum energy drift

0

Largest trajectory energy drift

0

Net outward window flux

0

Quantum / mean force / hopping milliseconds

0

Net outward flux subtracts incoming flux through the left window. Initially the packet is outside the windows: final exterior probability plus absorbed loss equals initial exterior probability plus net outward flux. Trigonometric interpolation evaluates current at the exact window planes; fractional-cell quadrature integrates their densities.

Costs measure propagation on this device, excluding snapshots, figures and downloads. Quantum energy drift is meaningful before absorption; the absorbing layer removes probability beyond nineteen bohr.

Model and numerical limits

The quantum solver uses 512 periodic Fourier points on a forty-eight-bohr box and a Strang time step no larger than 0.8 atomic time units. It applies the exact local two-state exponential and an outer absorbing mask. No stepwise normalization is applied. The offline benchmark halves the step and doubles the grid.

Ehrenfest trajectories follow the electronic mean force. Fewest-switches trajectories follow one adiabatic surface, propagate electronic amplitudes coherently, rescale momentum at accepted hops and reject energetically forbidden hops without velocity reversal. No decoherence correction is added.

Read the benchmark method

Start with the essentials

Focus question: Can more trajectories remove an approximation’s error in a reaction probability?
One-sentence intuition: More trajectories reduce sampling scatter. A systematic difference from the quantum solution can remain.

Core mathematical model

A nuclear coordinate coupled to two electronic states

H=-\frac{1}{2M}\partial_x^2 I+\begin{pmatrix}z(x)&w(x)\\w(x)&-z(x)\end{pmatrix},\quad M=2000

Atomic units set the reduced Planck constant to one. Position is in bohr, energy in hartree and time in atomic time units. The matrix eigenvalues are the two adiabatic energies.

The original simple avoided crossing

z(x)=0.01\,\operatorname{sgn}(x)(1-e^{-1.6|x|}),\qquad w(x)=0.005c\,e^{-x^2}

Electronic mixing multiplies the original off-diagonal coupling. The original Tully model has mixing equal to one. The initial packet is on the lower adiabatic surface.

A shared Gaussian preparation

\chi(x)\propto e^{-(x+8)^2/(4\sigma_x^2)+ip_0x},\quad \sigma_p=\frac{1}{2\sigma_x}

The classical methods independently sample the positive Wigner distribution of this same minimum-uncertainty packet. Their initial electronic state is the local lower eigenvector.

Common difficulties

A population is not an exit probability

Typical misconceptionElectronic population alone gives the reaction yield.

Better mental modelAn exit also requires the nuclear packet to reach its chosen spatial window. During the animation the incoming packet is still on the left; the table uses the common final time and reports unresolved weight separately.

Run the experiment

01
Watch the split
Choose Two outgoing branches and play the evolution. Compare the three rows at the same time.
What to observe: The quantum density spreads onto both energy curves. Surface-hopping trajectories separate by their active surface; mean-force trajectories can travel between the two curves.
02
Close an exit
Choose An upper exit closes and inspect the late-time channel probabilities.
What to observe: The upper asymptote requires enough total energy. A finite-width quantum packet contains an energy distribution; it can leave a small high-energy tail.
03
Separate scatter and bias
Open Ensemble and seed, increase Trajectories and save several independent seeds. Compare the final estimates with their standard errors.
What to observe: Repeated estimates fluctuate less when more trajectories are used. The gap between each method and the quantum result need not shrink.

Benchmark Nonadiabatic Reaction Methods

Separate sampling scatter from method bias

Four exit probabilities

Saved estimates and their sampling error