Part 0 introduced MQC as a closure problem: the electrons have amplitudes, but the nuclei need a force. Fewest-switches surface hopping keeps each trajectory on one active adiabatic surface, then uses stochastic hops to make the ensemble follow electronic population transfer with as few switches as possible.

Active-Surface Picture

For one trajectory, the electronic state is still propagated in the adiabatic basis,

\[ \dot c_j= -\frac{i}{\hbar}E_j(R)c_j -\sum_k \dot R\cdot d_{jk}(R)c_k. \]

The nuclei, however, feel only the active state \(a\):

\[ \dot R=\frac{P}{M}, \qquad \dot P=-\nabla_R E_a(R). \]

Electronic amplitudes and active surfaces are therefore different bookkeeping layers. The coefficients describe coherent population flow; the active state determines the classical force used by that trajectory.

Fewest-Switches Rule

Let \(\rho_{ij}=c_i c_j^*\). Differentiate the active-state population using the coefficient equation above. The diagonal energy term only changes phase and cancels from the population derivative:

\[ \begin{aligned} \dot\rho_{aa}&=2\operatorname{Re}(c_a^*\dot c_a)\\ &=-\sum_{b\ne a}2\operatorname{Re} (c_a^*c_b\,\dot R\cdot d_{ab}). \end{aligned} \]

Under the real-basis convention used here, the real part is unchanged by conjugating the electronic coherence. Thus the outgoing population flux can be written in the script's notation as

\[ J_{a\to b}=2\,\mathrm{Re}\!\left[ \rho_{ab}\,\dot R\cdot d_{ab} \right]. \]

If a fraction \(\rho_{aa}\) of otherwise equivalent trajectories occupies state \(a\), switching a fraction \(g_{a\to b}\) transfers \(\rho_{aa}g_{a\to b}\). Matching this to \(J_{a\to b}\Delta t\), and retaining only outward flow, gives the attempted hop probability:

\[ g_{a\to b}= \max\left(0,\frac{\Delta t\,J_{a\to b}}{\rho_{aa}}\right). \]

If a hop is selected, the nuclear momentum is rescaled so that the classical total energy is conserved:

\[ \frac{P'^2}{2M}+E_b(R) = \frac{P^2}{2M}+E_a(R). \]

For a one-dimensional accepted hop this gives \(P'=\operatorname{sgn}(P)\sqrt{P^2-2M(E_b-E_a)}\). If the kinetic energy is insufficient for an upward hop, the hop is frustrated and the trajectory remains on its original surface. In multidimensional simulations, this rescaling is normally done along the nonadiabatic-coupling direction.

Tully SAC Model

The test uses Tully's simple avoided crossing in a diabatic basis,

\[ V(x)= \begin{pmatrix} V_{11}(x) & V_{12}(x)\\ V_{12}(x) & -V_{11}(x) \end{pmatrix}, \]
\[ V_{11}(x)= \begin{cases} A[1-\exp(-Bx)], & x>0,\\ -A[1-\exp(Bx)], & x<0, \end{cases} \qquad V_{12}(x)=C\exp(-Dx^2), \]

with \(A=0.01\), \(B=1.6\), \(C=0.005\), \(D=1.0\), \(M=2000\), and \(\hbar=1\). The trajectory starts on the lower adiabatic state at \(x=-10\). The benchmark scans initial momentum \(p_0=10,14,18,20,22,26,30\), runs 50 trajectories at each momentum, and stops once the trajectory leaves the interaction region.

DVR Reference

The black reference curve is not another MQC trajectory. It is obtained by diagonalizing a two-channel sinc-DVR Hamiltonian on a finite grid,

\[ H_{\mathrm{DVR}}= \begin{pmatrix} T_{\mathrm{DVR}}+V_{11}(x_i)\delta_{ij} & V_{12}(x_i)\delta_{ij}\\ V_{12}(x_i)\delta_{ij} & T_{\mathrm{DVR}}+V_{22}(x_i)\delta_{ij} \end{pmatrix}, \]

then propagating an initial Gaussian wavepacket centered at \(x_0=-10\), width \(\sigma=1\), and mean momentum \(p_0\) by exact phase factors in that finite Hilbert space. The final wavepacket is projected onto local adiabatic states and integrated on the transmitted and reflected sides. This gives an exact quantum reference for the chosen DVR grid and initial packet, not a point-particle limit.

Result

In this momentum window every FSSH trajectory exits on the transmitted side, so the useful diagnostic is not reflection but which adiabatic branch is occupied after the crossing. At \(p_0=20\), the final FSSH branch probabilities are \(T_0=0.560\) and \(T_1=0.440\), while the DVR reference gives \(T_1=0.493\) with negligible reflection. As \(p_0\) increases, the upper transmitted branch becomes more likely; the mean final electronic upper population rises from 0.166 at \(p_0=10\) to 0.717 at \(p_0=30\). The residual oscillation in the active-state probabilities is finite-ensemble noise from the compact 50-trajectory run.

Fewest-switches surface hopping benchmark on Tully simple avoided crossing
FSSH benchmark on Tully's simple avoided crossing. Left: horizontal axis is nuclear coordinate \(x\); vertical axis is energy in atomic units, with the derivative coupling scaled onto the same axis. Blue and red dashed curves are diabatic \(V_{11}\) and \(V_{22}\), black and gray solid curves are adiabatic energies \(E_0\) and \(E_1\), and the green curve is \(d_{01}/160\). Right: horizontal axis is initial momentum \(p_0\); vertical axis is probability. Blue circles are FSSH transmitted lower-surface active-state probability, red squares are FSSH transmitted upper-surface probability, gray triangles are FSSH total reflection probability, black diamonds are the DVR exact transmitted upper probability for the Gaussian packet, gray inverted triangles are the DVR exact reflection probability, and the green dotted curve is the mean final electronic upper population \(|c_1|^2\). The main comparison is that FSSH tracks the DVR upper-branch trend while retaining stochastic active-state noise.

Code used in this note

  • fssh_tully_sac.py runs the FSSH ensemble, applies energy-conserving hops, writes CSV tables, and renders the benchmark figure.
  • dvr_tully_sac_reference.py builds and diagonalizes the finite sinc-DVR two-channel Hamiltonian used for the exact quantum reference.
  • tully_common.py defines the Tully SAC potential, adiabatic energies, forces, derivative coupling, and shared electronic propagation helpers.
  • fssh-tully-sac.csv stores the momentum scan, branch probabilities, mean hop counts, and final electronic populations.
  • fssh-tully-sac-summary.csv stores the compact numerical values quoted in the text.
  • dvr-tully-sac-reference.csv stores the DVR transmitted/reflected lower/upper adiabatic probabilities.

References

  • Tully, J. Chem. Phys. 93, 1061 (1990).
  • Crespo-Otero and Barbatti, Chem. Rev. 118, 7026 (2018).
  • Nelson, White, Bjorgaard, Sifain, Zhang, Nebgen, Fernandez-Alberti, Mozyrsky, Roitberg, and Tretiak, Chem. Rev. 120, 2215 (2020).

MQC series map

  1. Part 0: Background and Working Assumptions - BO approximation, Born-Huang couplings, NAC/HST formulas, MQC closures, and limitations.
  2. Part I: Fewest-Switches Surface Hopping - active surfaces, stochastic hops, velocity rescaling, and a Tully SAC branch benchmark.
  3. Part II: Ehrenfest Mean-Field Dynamics - mean-field forces, continuous electronic populations, and the same Tully SAC scan.
  4. Part III: P-Matrix Decoherence and Detailed Balance - directed density-matrix variables, pairwise decoherence, Boltzmann detailed balance, and controlled toy diagnostics.
  5. Part IV: MMST Mapping and ZPE Leakage - oscillator mapping variables, electronic phase sampling, ZPE leakage, gamma choices, SQC-style corrections, and spin-mapping context.