Mixed quantum-classical dynamics is a practical compromise for nonadiabatic molecular motion: electrons are propagated quantum mechanically, while nuclei are represented by classical trajectories. This opening note sets up the common equations and vocabulary used in the following FSSH and Ehrenfest benchmarks. It is intentionally conceptual; executable tests start in Parts I and II.

Why MQC Is Needed

The Born-Oppenheimer picture assumes that the electronic state adapts instantaneously to nuclear motion. That separation becomes fragile near avoided crossings, conical intersections, charge-transfer regions, spin-orbit crossings, and strong laser-field couplings. In those regions, nuclear motion can move population between electronic states, and a single potential energy surface is no longer enough.

Fully quantum nuclear-electronic dynamics is usually too expensive for realistic molecular systems. NA-MQC methods reduce the cost by keeping the most phase-sensitive subsystem, usually the electrons, in a quantum representation while replacing nuclear wavepackets by an ensemble of classical paths. This is why the methods are common in photochemistry, exciton migration, organic semiconductors, charge separation and recombination, radiation chemistry, and ultrafast spectroscopy.

Born-Oppenheimer Starting Point

Start from a molecular Hamiltonian split into nuclear kinetic energy and an electronic Hamiltonian parametrized by nuclear coordinates:

\[ \hat H=\hat T_n+\hat H_e(r;R), \qquad \hat H_e(r;R)\phi_j(r;R)=E_j(R)\phi_j(r;R). \]

The Born-Oppenheimer approximation assumes that a nuclear wavepacket can remain on one electronic eigenstate,

\[ \Psi(r,R,t)\approx \chi_\alpha(R,t)\phi_\alpha(r;R), \qquad i\hbar\frac{\partial\chi_\alpha}{\partial t} = \left[\hat T_n+E_\alpha(R)\right]\chi_\alpha . \]

This approximation is useful because the nuclei move on a single potential energy surface. It fails when the electronic basis changes rapidly along the nuclear path, especially near small energy gaps. MQC methods can be viewed as controlled attempts to keep some of the dropped coupling while avoiding fully quantum nuclear dynamics.

Born-Huang Couplings

The exact adiabatic expansion is

\[ \Psi(r,R,t)=\sum_j \chi_j(R,t)\phi_j(r;R). \]

The nuclear kinetic operator differentiates both the nuclear coefficient and the electronic basis. For one coordinate \(R\), the product rule gives

\[ \partial_R^2(\chi_j\phi_j) =\phi_j\chi_j''+2\phi_j'\chi_j'+\phi_j''\chi_j. \]

Insert the expansion into the molecular Schrödinger equation, multiply by \(\phi_i^*\), and integrate over electronic coordinates. Orthonormality selects the ordinary kinetic term, while the other two derivatives produce the coupled-channel structure:

\[ i\hbar\dot\chi_i = \left[\hat T_n+E_i(R)\right]\chi_i -\sum_j\frac{\hbar^2}{2M} \left[ 2d_{ij}(R)\frac{\partial}{\partial R} +D_{ij}(R) \right]\chi_j, \]
\[ d_{ij}(R)=\langle\phi_i|\nabla_R\phi_j\rangle, \qquad D_{ij}(R)=\langle\phi_i|\nabla_R^2\phi_j\rangle . \]

The strict BO approximation drops the off-diagonal channel couplings and often absorbs or neglects the diagonal correction. Nonadiabatic dynamics begins when those derivative terms cannot be ignored.

NAC Working Formulas

The first-derivative coupling, or NAC vector, is

\[ d_{ij}(R)=\langle\phi_i(r;R)|\nabla_R\phi_j(r;R)\rangle. \]

To obtain a working expression, differentiate \(\hat H_e\phi_j=E_j\phi_j\) and project onto a different eigenstate. The derivative of \(E_j\) drops out by orthogonality:

\[ \langle\phi_i|\nabla_R\hat H_e|\phi_j\rangle +E_i d_{ij}=E_jd_{ij},\qquad i\ne j. \]

For nondegenerate adiabatic states, rearrangement gives the Hellmann–Feynman form:

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

When the energy gap becomes small, derivative coupling can become large. This is the mathematical signal that the nuclei can no longer be assumed to remain on one isolated adiabatic surface.

Along one classical trajectory, codes often need the scalar time-derivative coupling

\[ \tau_{ij}(t)=\langle\phi_i|\dot\phi_j\rangle = \dot R(t)\cdot d_{ij}(R(t)). \]

If analytic NAC vectors are unavailable, the Hammes-Schiffer-Tully overlap formula estimates this quantity from electronic wavefunctions at two adjacent nuclear geometries. With the overlap matrix

\[ S_{ij}(t,t+\Delta t)= \langle\phi_i(t)|\phi_j(t+\Delta t)\rangle, \]

the antisymmetric finite-difference estimate is

\[ \tau_{ij}\!\left(t+\frac{\Delta t}{2}\right) \approx \frac{S_{ij}(t,t+\Delta t)-S_{ji}^*(t,t+\Delta t)} {2\Delta t}, \]

which reduces to the usual antisymmetric real-state expression after phase matching. In complex or spin-orbit-coupled calculations, the conjugated overlap convention must be handled carefully. The formula is practical, but it is not a substitute for state tracking: arbitrary sign flips, root reordering, and near-degenerate subspaces can turn a clean NAC estimate into numerical noise.

Classical Nuclei

MQC replaces the nuclear wavepacket by trajectories \((R(t),P(t))\), while the electronic state is expanded along each moving path:

\[ |\psi_e(t)\rangle=\sum_j c_j(t)|\phi_j(R(t))\rangle. \]

Differentiate both the coefficients and the moving basis. Since \(\dot\phi_k=\dot R\cdot\nabla_R\phi_k\), projection onto state \(j\) gives

\[ i\hbar\left(\dot c_j+ \sum_k\dot R\cdot d_{jk}\,c_k\right)=E_jc_j. \]

Rearranging gives the coefficient equation used by the trajectory methods:

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

The first term is phase accumulation on the adiabatic energy; the second term transfers amplitude when the trajectory crosses a region with nonzero derivative coupling. The remaining question is how the classical nuclei feel the electronic state. Different MQC methods answer that closure differently.

Two Closures

Ehrenfest dynamics uses a mean force. In a two-state adiabatic model, the force has the schematic form

\[ F_{\mathrm{Eh}} = \rho_{00}F_0+\rho_{11}F_1 +2\,\mathrm{Re}\!\left[\rho_{01}d_{01}\right](E_1-E_0), \qquad \rho_{ij}=c_i c_j^*. \]

This is smooth and deterministic, but a single trajectory cannot split into distinct nuclear branches after the electronic wavepacket separates.

Fewest-switches surface hopping instead keeps one active adiabatic surface for each trajectory. The electronic amplitudes still evolve continuously, but nuclear forces come from the active state. Stochastic hops are chosen so that the ensemble of active states follows the electronic population flow as closely as possible. This gives a branching picture, but it introduces practical questions: velocity rescaling, frustrated hops, decoherence corrections, detailed balance, and ensemble convergence.

Useful Domain

MQC is most useful when electronic transitions are essential but the nuclear motion is sufficiently classical over the time scale of interest. It is a working language for direct dynamics because forces, energies, and nonadiabatic couplings can be evaluated locally along trajectories rather than precomputing global potential energy surfaces.

The limitations are equally important. Classical nuclei miss nuclear tunneling, zero-point leakage control, nuclear interference, and exact wavepacket splitting. Electronic coherence often persists too long unless a decoherence correction or a more elaborate trajectory-coupled method is used. The following notes use Tully's one-dimensional avoided-crossing models precisely because they expose these issues without the clutter of a full electronic-structure calculation.

References

  • 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).
  • Tully, J. Chem. Phys. 93, 1061 (1990).
  • Hammes-Schiffer and Tully, J. Chem. Phys. 101, 4657 (1994).

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.