After the DVR Hamiltonian is built, the next step is a real-time wavepacket calculation. The initial state is a Gaussian nuclear packet placed on one electronic state, then propagated by the full product-basis Hamiltonian.

What the propagation is

The calculation is not an iterative DVR solver for the Schrodinger equation. DVR first turns the continuous coordinate problem into one finite Hamiltonian matrix. For the time-independent models here, the evolution at any requested time is then evaluated directly from the spectral decomposition of that matrix:

Build the product-basis Hamiltonian
\(H_{(i,a),(j,b)}=\delta_{ab}T_{ij}+\delta_{ij}V_{ab}(x_i)\)
Diagonalize once
\(H=C\,\mathrm{diag}(E_\alpha)C^\dagger\)
Evaluate an arbitrary saved time
\(\psi(t)=C\,\mathrm{diag}(e^{-iE_\alpha t/\hbar})C^\dagger\psi(0)\)
Read observables
Sum \(|\psi_a(x_i,t)|^2\Delta x\) for populations, densities, and packet diagnostics.

Loops in the plotting code are therefore loops over requested output times. They are not Euler, Runge-Kutta, or finite-difference time steps approximating the time derivative.

Product basis

For \(N_x\) DVR grid points and \(n_{\mathrm{state}}\) electronic states, the working basis is

\[ |x_i,a\rangle,\qquad i=1,\ldots,N_x,\quad a=1,\ldots,n_{\mathrm{state}}. \]

The Hamiltonian has the same block form as in the previous note:

\[ H_{(i,a),(j,b)}=\delta_{ab}T_{ij}+\delta_{ij}V_{ab}(x_i). \]

The diabatic potential \(V(x_i)\) is evaluated independently at every grid point. Local adiabatic surfaces are obtained by diagonalizing this \(n_{\mathrm{state}}\times n_{\mathrm{state}}\) matrix at each \(x_i\).

Initial Gaussian packet

The nuclear part of the initial wavefunction is

\[ \chi(x,0)=A \exp\!\left[ik_0(x-x_0)\right] \exp\!\left[-\sigma(x-x_0)^2\right], \]

where \(x_0\) fixes the packet center, \(k_0\) fixes the mean momentum, and \(\sigma\) controls the spatial width. On the DVR grid it becomes a vector, normalized by

\[ \sum_i |\chi(x_i,0)|^2\Delta x=1. \]

If the packet is prepared on the lower local adiabatic state, the source notes transform it into the diabatic product basis before propagation.

Potential model and state populations

Tully extended coupling reflection potential
Example two-state potential. The blue curve is \(V_{11}(x)\), the orange curve is \(V_{22}(x)\), and the green dashed curve is the diabatic coupling \(V_{12}(x)\). The black and gray solid curves are the local adiabatic surfaces \(E_1(x)\) and \(E_2(x)\), read on the left energy axis. The red dotted curve is the nonadiabatic derivative coupling magnitude \(|d_{12}(x)|\), read on the right axis; its peak marks the region where state transfer is strongest.
Population dynamics from DVR wavepacket propagation
Population dynamics from the propagation. The blue curve is state 1 population and the orange curve is state 2 population. The horizontal axis is time in atomic units; the vertical axis is total probability on each electronic state.

Propagation by spectral decomposition

For a time-independent Hamiltonian, the finite-dimensional evolution is

\[ \psi(t)=U(t)\psi(0),\qquad U(t)=e^{-iHt/\hbar}. \]

After diagonalizing \(H=C\,\mathrm{diag}(E_\alpha)C^\dagger\), this becomes

\[ \psi(t)= C\,\mathrm{diag}\!\left(e^{-iE_\alpha t/\hbar}\right)C^\dagger\psi(0). \]

Electronic-state populations are then simple sums over grid points:

\[ P_a(t)=\sum_i |\psi_a(x_i,t)|^2\Delta x. \]

Within the chosen finite basis, this spectral propagation is exact up to diagonalization and floating-point roundoff. The approximation enters earlier: in the finite coordinate basis, the box size, the grid spacing, and the potential model.

Density heatmaps

Nuclear density heatmap from wavepacket propagation
Nuclear density heatmap. The horizontal axis is nuclear coordinate \(x\); the vertical axis is saved time \(t\). Color represents \(\rho(x,t)=|\psi(x,t)|^2\), with dark purple near zero and yellow at high density. Bright bands trace where the packet is concentrated as it moves through the coupling region.
Long-time nuclear density heatmap showing boundary effects
Longer propagation diagnostic. The horizontal axis is coordinate and the vertical axis is time. The color scale is \(\log_{10}\rho(x,t)\), so yellow indicates larger density and purple indicates very small density. The blue overlaid line tracks the mean position \(\langle x\rangle(t)\). Its turn near the finite-box edge shows why boundary handling becomes important in later calculations.

From wavepackets to ensembles

The single-packet calculation prepares a pure state. Part IV keeps the same finite Hamiltonian idea but replaces \(\psi(t)\) with a density matrix \(\rho(t)\), which makes thermal ensembles and incoming momentum filters natural.

Code used in this note

  • dvr_kubo_minimal.py includes the two-state Tully Hamiltonian constructor, Gaussian packet initializer, and spectral propagation routine used for this note.
  • source/README.md describes the cleaned code attachments and what was kept out of the published site.

DVR series map

  1. Part I: Wavepacket Dynamics Methods - discretization choices and the harmonic-oscillator benchmark.
  2. Part II: Operator Matrices in the DVR Basis - matrix structure behind kinetic, side, and flux operators.
  3. Part III: Wavepacket Dynamics with DVR - exact finite-dimensional propagation of a prepared wavepacket.
  4. Part IV: Density-Matrix Ensemble Evolution - thermal ensembles, filtering, and absorber bookkeeping.
  5. Part V: Kubo-Transformed Position Correlation - spectral thermal correlation functions.
  6. Part VI: Flux-Side Kubo Correlation - rate-like observables from side and flux operators.
  7. Part VII: Excited-Ground Kubo Population Correlation - local adiabatic projectors and population-resolved transfer diagnostics.