Fewest Switches Surface Hopping

• FSSH is a quantum–classical algorithm that simulates nonadiabatic dynamics by combining classical trajectory propagation with stochastic surface hops.
• It employs time-dependent electronic propagation and energy-conserving momentum rescaling to accurately reproduce quantum state populations in complex systems.
• Ongoing extensions incorporate decoherence corrections, multistate and time-dependent enhancements, and machine learning approaches to improve simulation accuracy.

Fewest Switches Surface Hopping (FSSH) is the dominant mixed quantum–classical algorithm for simulating nonadiabatic dynamics in complex molecular and condensed-phase systems. It combines classical trajectory propagation on adiabatic potentials with quantum evolution of electronic amplitudes and implements stochastic “hops” between surfaces to reproduce quantum populations with minimal switching. FSSH is widely used for simulating radiationless decay, photochemical reactions, energy transfer, and electronic relaxation in molecules and materials, with applications extending from solution-phase dynamics to nanostructures and crystalline solids.

1. Formal Foundation: QCLE Perspective and Algorithmic Structure

FSSH is grounded in the quantum–classical Liouville equation (QCLE), which rigorously defines the evolution of the partially Wigner-transformed electronic density matrix $\rho_W^{\alpha\alpha'}(X,t)$ as a function of nuclear coordinates and momenta $X=(R,P)$. In the adiabatic basis, the QCLE reads

$$\frac{\partial}{\partial t}\,\rho_W^{\alpha\alpha'}(X,t) = -\,i\omega_{\alpha\alpha'}(R)\,\rho_W^{\alpha\alpha'}(X,t) -\;i\,L_{\alpha\alpha'}\,\rho_W^{\alpha\alpha'}(X,t) +\sum_{\beta\beta'}\mathcal J_{\alpha\alpha',\beta\beta'}(X)\,\rho_W^{\beta\beta'}(X,t),$$

where $$\omega_{\alpha\alpha'}(R) = [E_\alpha(R) - E_{\alpha'}(R)]/\hbar$$, the classical Liouville operator $L_{\alpha\alpha'}$ propagates $X$ on either the ground or the mean surface, and the nonadiabatic coupling superoperator $\mathcal J$ encodes both population transfer (nonadiabatic hops) and quantum “kicks” (momentum changes) (Kapral, 2016).

FSSH emerges by applying three key approximations:

1. Neglect explicit decoherence by setting the QCLE decoherence operator $\gamma_{\alpha\alpha'}^{(\nu)}$ to zero.
2. Restrict momentum jumps to transitions between populations (diagonal density matrix elements), fusing infinitesimal QCLE kicks into finite, energy-conserving momentum rescalings.
3. Drop trajectory reweighting—hopping probabilities are chosen to give positive-definite populations without accounting for quantum branching weights.

The resultant FSSH algorithm propagates an ensemble of trajectories, each initialized on a single adiabatic surface with classical nuclei. At each step, the electronic amplitudes $c_j$ are evolved using the time-dependent Schrödinger equation,

$$i\hbar\,\dot{c}_j = E_j(R)\,c_j - i\hbar\sum_k (P/M)\cdot d_{jk}(R)\,c_k,$$

and stochastic hops are introduced between surfaces, with probabilities designed to match the ensemble-averaged adiabatic populations to $|c_i|^2$ as closely as possible (Kapral, 2016, Bondarenko et al., 2022).

2. Electronic Propagation, Hopping Probabilities, and Momentum Rescaling

Each trajectory carries a vector of complex electronic amplitudes $\{c_k(t)\}$ over adiabatic states $\{|k\rangle\}$ and evolves them according to the time-dependent Schrödinger equation in the instantaneous basis. The hopping probability for transition $k \rightarrow l$ during a time interval $[t, t+\Delta t]$ is

$$g_{k\to l} = 2\,\mathrm{Re}\left[\dot{q}\cdot d_{kl} \;\frac{c_l\,c_k^*}{|c_k|^2}\right]\,\Delta t,$$

subject to $0\leq g_{k\to l}\leq 1$ per step (Bondarenko et al., 2022). Hops are enacted stochastically using random numbers, enforcing the “fewest switches” principle.

Upon a successful hop, classical momentum is rescaled to conserve total energy. The canonical FSSH momentum rescaling along the unit nonadiabatic coupling direction $\hat d_{kl}$ is

$$\Delta P_{kl}^{FS} = \hat d_{kl} \, \mathrm{sgn}(P \cdot \hat d_{kl}) \sqrt{(P \cdot \hat d_{kl})^2 + 2M (E_k - E_l)} - P \cdot \hat d_{kl},$$

with the hop rejected if the initial kinetic energy along $\hat d_{kl}$ does not suffice for energy conservation. Recent first-principles derivations from quantum trajectory surface hopping (QTSH) confirm this formula in the impulsive, localized-transition limit, linking FSSH rescaling to the integrated quantum force that transfers population between adiabats (Huang et al., 2023).

3. Decoherence, Detailed Balance, and Algorithmic Variations

A central limitation of unconstrained FSSH is the lack of dynamical decoherence: electronic coherences remain finite even as nuclear wave packets on different surfaces diverge, leading to overestimated coherence lifetimes, incorrect quantum branching ratios, and violations of detailed balance in certain regimes (Kapral, 2016, Chen et al., 2016, Nagda et al., 22 Nov 2025). The QCLE formalism provides an explicit decoherence rate,

$$\gamma_{\alpha\alpha'}^{(\nu)}(X) = \frac{1}{2}\left[F_\alpha(R) - F_\nu(R) + F_{\alpha'}(R) - F_\nu(R)\right] \cdot \frac{\partial_P \rho_W^{\alpha\alpha'}(X)}{\rho_W^{\alpha\alpha'}(X)}.$$

Practical decoherence corrections include:

• Augmented FSSH (A-FSSH): Estimates decoherence times based on wavepacket separation and collapses off-diagonal density matrix elements once the nuclear overlap vanishes, restoring correct scaling in biased golden rule rates and damping overcoherent oscillations (Chen et al., 2016, Nagda et al., 22 Nov 2025).
• Energy-based, Gaussian-overlap, or instantaneous collapse schemes: Apply local decoherence rates only in weak-coupling regions, exploiting the dimensionless Massey parameter $r = \hbar |T_{ab}|/|V_a - V_b|$ to switch decoherence on/off (Runeson, 24 Jul 2025).

Despite these advances, some regimes remain problematic. For example, A-FSSH yields correct rate constants but produces incorrect thermal populations (with an effective $T_{\mathrm{eff}}\approx T/2$) in the deep Marcus inverted regime due to self-consistency issues between trajectories and quantum populations (Nagda et al., 22 Nov 2025). This highlights the nuanced interplay between decoherence, detailed balance, and observable accuracy.

4. Extensions: Multi-State, Open and Driven Systems, Machine Learning, and Alternative Representations

The FSSH protocol has been adapted to various complex systems and regimes:

• Multistate and band-like systems: Real- and reciprocal-space formulations have been developed, extending surface hopping to crystalline or periodic systems with large band manifolds, and enabling direct integration with band-structure methods (Krotz et al., 2021).
• Driven quantum-classical dynamics: FSSH has been extended for explicit time-dependent Hamiltonians (e.g., laser-driven photodynamics), utilizing Floquet representations and quasienergy surfaces, and accounting for time-dependent nonadiabatic couplings (Heindl et al., 2021, Zhou et al., 2019).
• Complex Hamiltonians and geometric phases: FSSH has been generalized for systems with Berry curvature or magnetic fields, introducing gauge-invariant momentum rescaling, magnetic (Berry) forces, and explicit accounting of topological effects (Miao et al., 2019, Wu et al., 2023).
• Machine learned dynamics: Recent progress enables ML-prediction of energies, gradients, and nonadiabatic couplings via descriptors tailored for phase consistency and singularity removal, facilitating fully ML-driven FSSH with accurate reproduction of nonadiabatic population dynamics at massively reduced computational cost (Martinka et al., 29 May 2025, Tang et al., 2022).
• Numerical stability and time step adaptation: FSSH-2 eliminates explicit nonadiabatic coupling vectors, using local-diabatization overlap propagation and population-based hop probabilities, leading to robust dynamics even with large integration steps (Araujo et al., 18 Jan 2024).

5. Reduced Density Matrix Construction, Positivity, and Observables

Standard FSSH population estimates are inherently consistent with the fraction of trajectories on each state, but construction of a physically meaningful reduced density matrix—especially for observables dependent on coherences—poses challenges. Naïve combination of populations (from stochastic hops) and coherences (from wavefunction amplitudes) can violate positivity and lead to negative eigenvalues or unphysical observables, as exemplified in the homogeneous trimer model at low temperatures (Bondarenko et al., 2022).

Coherent FSSH (C-FSSH) constructs the density matrix exclusively from active-surface branch pairs, weighting off-diagonal elements with phase accumulations and classical overlaps,

$$\rho_{kl}^{(i,j)}(t) = \rho_0^{(i,j)}\,\delta_{k,a^{(i)}}\,\delta_{l,a^{(j)}}\,F^{(i,j)}(t)\,e^{-i[\phi_k^{(i)}-\phi_l^{(j)}]},$$

guaranteeing $\rho\ge0$ for all parameter regimes tested. C-FSSH typically improves population accuracy and preserves positivity, though can overdamp coherences in the Marcus regime (Bondarenko et al., 2022).

6. Practical Implementation, Benchmarking, and Limitations

FSSH has been implemented in a wide array of electronic structure–dynamics frameworks and ab initio codes. The workflow typically involves:

• Generating initial conditions via Wigner or Boltzmann sampling,
• On-the-fly computation of energies, forces, and nonadiabatic couplings,
• Propagation of nuclear and electronic DOFs,
• Stochastic surface transitions and momentum adjustments,
• Averaging over large trajectory ensembles for statistical observables (Miller et al., 15 Feb 2024).

Extensive benchmarking reveals:

• FSSH generally yields accurate nonadiabatic population and rate predictions across weak–intermediate coupling, high-temperature, and large-system limits (Chen et al., 2016, Krotz et al., 2021).
• Decoherence corrections are essential for detailed balance in biased regimes, quantum yields, and correct suppression of spurious wavepacket recurrences (Chen et al., 2016, Runeson, 24 Jul 2025).
• Explicit DBOC inclusion is typically detrimental near conical intersections or for strong couplings, and should be omitted except in weak-coupling, low-kinetic-energy regimes (Gherib et al., 2016).
• For spin–orbit coupled and chiral systems, phase-space FSSH methods (PSSH) constructed on phase-space Hamiltonians $H(R,P)$ with electronic–nuclear couplings $\Gamma\cdot P$ can recover exact conservation of total linear and angular momentum, essential for simulating chiral-induced spin selectivity (CISS) and related effects (Wu et al., 2023).

A summary of main FSSH properties and challenges is provided below.

Aspect	Standard FSSH	Notable Limitations / Remedies
Population Accuracy	Good in most regimes	Overcoherence errors, especially at long times
Coherence/Decoherence	Lacks dynamical decoherence	A-FSSH, GONT, INT, C-FSSH schemes
Density Matrix	Positivity violation common	C-FSSH: strictly positive-definite $\rho$
Band/Extended Systems	Real-space FSSH inefficient	Reciprocal-space FSSH for band structure
Machine Learning	Typically QM/ML or hybrid implementations	Phase correction, tailored descriptors required
DBOC inclusion	Not present	Detrimental near CIs or strong coupling regimes

7. Outlook and Future Directions

The FSSH framework continues to evolve, with advances targeting:

• Rigorous incorporation of decoherence and thermalization,
• Parameter-free density matrix construction,
• Efficient ML-based surrogates for quantum observables and nonadiabatic couplings,
• Extensions for time-dependent, magnetic, and chiral Hamiltonians,
• Adaptive and energy-conserving algorithms for systems with geometric and topological structure.

Ongoing research aims to bridge the gap to exact quantum-classical and fully quantum benchmarks, increase computational tractability for extended systems, and refine the interpretation of nonequilibrium and time-resolved experimental observables. The FSSH paradigm remains a central tool for predictive nonadiabatic molecular and materials simulations (Kapral, 2016, Bondarenko et al., 2022, Huang et al., 2023).