.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/semiclassical/propagators/plot_heller_franck_condon_spectroscopy.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_semiclassical_propagators_plot_heller_franck_condon_spectroscopy.py: Heller's time-dependent spectroscopy: absorption bands from short-time dynamics ================================================================================== Heller (1981) recast molecular absorption spectroscopy in the time domain. A photon lifts the vibrational ground state :math:`\psi_0` of the lower electronic surface onto the upper surface, where it is no longer stationary. The absorption spectrum is the Fourier transform of the autocorrelation of the packet that then moves on the upper surface, .. math:: \sigma(E) \propto \int_{-\infty}^{\infty}\langle\psi_0|\psi(t)\rangle\, e^{iEt/\hbar}\,dt = 2\pi\hbar\sum_n|\langle n|\psi_0\rangle|^2\,\delta(E-E_n), the Franck-Condon factors of the upper surface's levels. The time-domain view explains the spectrum by events of the packet's motion. Its departure from the Franck-Condon region, within a fraction of a period, sets the width of the whole band. Its return after a vibrational period, if the upper surface is bound, splits the band into a progression. The individual lines need the long-time signal. The first two stages take a period or two, so short-time semiclassical dynamics already accounts for the band shape and the progression. This example uses one thawed Gaussian (:func:`~physicskit.semiclassical.core.propagators.thawed_gaussian_propagate`) on two upper surfaces, a bound Morse curve and a repulsive one, and compares the spectra with the exact Franck-Condon factors (:class:`~physicskit.quantum.core.eigensolvers.NumerovSolver`) and with split-operator evolution (:class:`~physicskit.quantum.core.solvers.SplitOperatorSolver1D`). Units: :math:`\hbar=m=1`. .. GENERATED FROM PYTHON SOURCE LINES 36-119 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from numba import njit from physicskit.quantum.core.eigensolvers import NumerovSolver from physicskit.quantum.core.solvers import SplitOperatorSolver1D from physicskit.semiclassical.core.propagators import frozen_gaussian_1d, thawed_gaussian_propagate, thawed_gaussian_wavefunction # Ground surface: harmonic, omega = 2, so the initial packet is a Gaussian # with gamma = omega / 2 = 1 centred at x = 0. gamma0 = 1.0 x = np.linspace(-8, 88, 8192) dx = x[1] - x[0] psi0 = frozen_gaussian_1d(x, 0.0, 0.0, gamma0) # Bound upper surface: Morse well displaced to x = d, raised by T_e. D_e, a_m, d, T_e = 40.0, 0.25, 2.0, 10.0 morse_params = np.array([D_e, a_m, d, T_e]) @njit def V_morse(q, params): return params[3] + params[0] * (1 - np.exp(-params[1] * (q - params[2]))) ** 2 @njit def dV_morse(q, params): e = np.exp(-params[1] * (q - params[2])) return 2 * params[0] * params[1] * e * (1 - e) @njit def d2V_morse(q, params): e = np.exp(-params[1] * (q - params[2])) return 2 * params[0] * params[1] ** 2 * e * (2 * e - 1) # Repulsive upper surface: an exponential wall. A_r, b_r = 16.0, 0.6 rep_params = np.array([A_r, b_r]) @njit def V_rep(q, params): return params[0] * np.exp(-params[1] * q) @njit def dV_rep(q, params): return -params[0] * params[1] * np.exp(-params[1] * q) @njit def d2V_rep(q, params): return params[0] * params[1] ** 2 * np.exp(-params[1] * q) omega_e = a_m * np.sqrt(2 * D_e) T_vib = 2 * np.pi / omega_e dt = 0.004 n_steps = int(round(6 * T_vib / dt)) save = 5 def autocorrelations(V, dV, d2V, params): _, q, p, alpha, s = thawed_gaussian_propagate(0.0, 0.0, gamma0, dV, d2V, V, 1.0, dt, n_steps, params=params) ks = np.arange(0, n_steps + 1, save) C_tg = np.array([np.sum(psi0.conj() * thawed_gaussian_wavefunction(x, q[k], p[k], alpha[k], s[k])) * dx for k in ks]) solver = SplitOperatorSolver1D(x, lambda xx: V(xx, params), dt=dt) frames, times = solver.propagate(psi0, n_steps, save_every=save) C_ex = np.array([np.sum(psi0.conj() * f) * dx for f in frames]) return times, C_ex, C_tg def spectrum(E, times, C, t_window): w = np.exp(-0.5 * (times / t_window) ** 2) dts = times[1] - times[0] return (2 * np.real(np.exp(1j * np.outer(E, times)) @ (C * w)) - (C[0] * w[0]).real) * dts / (2 * np.pi) times, C_morse_ex, C_morse_tg = autocorrelations(V_morse, dV_morse, d2V_morse, morse_params) _, C_rep_ex, C_rep_tg = autocorrelations(V_rep, dV_rep, d2V_rep, rep_params) .. GENERATED FROM PYTHON SOURCE LINES 120-125 The two upper surfaces -------------------------- The vertical transition puts the packet on the steep inner side of the displaced Morse curve, or on the repulsive wall. Either way it starts to slide downhill at once. .. GENERATED FROM PYTHON SOURCE LINES 125-138 .. code-block:: Python fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4)) xx = np.linspace(-2.5, 9, 400) ax1.plot(xx, 2 * xx**2, color="0.5", label="ground surface") ax1.plot(xx, [V_morse(v, morse_params) for v in xx], color="steelblue", label="bound upper surface (Morse)") ax1.plot(xx, [V_rep(v, rep_params) for v in xx], color="firebrick", label="repulsive upper surface") ax1.fill_between(xx, 1.0, 1.0 + 3 * np.abs(frozen_gaussian_1d(xx, 0, 0, gamma0)) ** 2, color="0.5", alpha=0.4) ax1.annotate("", (0, V_morse(0.0, morse_params)), (0, 1.1), arrowprops=dict(arrowstyle="->", color="k")) ax1.set_ylim(0, 30) ax1.set_xlabel("vibrational coordinate x") ax1.set_ylabel("energy") ax1.set_title("Vertical (Franck-Condon) excitation") ax1.legend(fontsize=8) .. image-sg:: /api/gallery/semiclassical/propagators/images/sphx_glr_plot_heller_franck_condon_spectroscopy_001.png :alt: Vertical (Franck-Condon) excitation :srcset: /api/gallery/semiclassical/propagators/images/sphx_glr_plot_heller_franck_condon_spectroscopy_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 139-148 Autocorrelation: departure and return ----------------------------------------- On both surfaces :math:`|C(t)|` first falls as the packet leaves, and the thawed Gaussian follows that decay closely. On the Morse surface the packet comes back a little more than one small-oscillation period later (the well softens with energy). The thawed Gaussian comes back at about the right time, but as a single, strongly stretched Gaussian it gets the shape of the return only roughly. On the repulsive surface the packet never returns. .. GENERATED FROM PYTHON SOURCE LINES 148-159 .. code-block:: Python ax2.plot(times / T_vib, np.abs(C_morse_ex), color="steelblue", lw=3, alpha=0.4, label="Morse, exact") ax2.plot(times / T_vib, np.abs(C_morse_tg), color="steelblue", lw=1, label="Morse, thawed Gaussian") ax2.plot(times / T_vib, np.abs(C_rep_ex), color="firebrick", lw=3, alpha=0.4, label="repulsive, exact") ax2.plot(times / T_vib, np.abs(C_rep_tg), color="firebrick", lw=1, label="repulsive, thawed Gaussian") ax2.set_xlim(0, 3) ax2.set_xlabel(r"t / $T_{\rm vib}$") ax2.set_ylabel(r"$|\langle\psi_0|\psi(t)\rangle|$") ax2.set_title("Autocorrelation on the upper surface") ax2.legend(fontsize=8) fig1.tight_layout() .. GENERATED FROM PYTHON SOURCE LINES 160-171 Spectra ----------- The exact Franck-Condon factors are drawn as sticks. A short time window, about a quarter period, sees only the departure and gives the smooth band. Heller's estimate of it is a Gaussian centred at the vertical energy :math:`V_e(0)` with width :math:`|V_e'(0)|\sigma_x`, the slope of the upper surface times the packet width. A window of about one period (:math:`\sigma_t=0.8\,T`) includes the first return and splits the band into the vibrational progression, which the thawed Gaussian reproduces near the band centre and less well towards its edges. For the repulsive surface there is no return and no progression. .. GENERATED FROM PYTHON SOURCE LINES 171-211 .. code-block:: Python E = np.linspace(0, 40, 2000) sigma_x = np.sqrt(1 / (4 * gamma0)) res = NumerovSolver(np.linspace(-4, 30, 3000), lambda r: V_morse(r, morse_params)).solve(n_states=16) fc = np.array([np.trapezoid(w * np.exp(-gamma0 * res.x**2) * (2 * gamma0 / np.pi) ** 0.25, res.x) ** 2 for w in res.wavefunctions]) fig2, (ax3, ax4) = plt.subplots(1, 2, figsize=(12, 4)) for t_win, color in [(0.25 * T_vib, "0.5"), (0.8 * T_vib, "steelblue")]: ax3.plot(E, spectrum(E, times, C_morse_ex, t_win), color=color, lw=3, alpha=0.4) ax3.plot(E, spectrum(E, times, C_morse_tg, t_win), color=color, lw=1, ls="--", label=f"thawed, window {t_win / T_vib:.2f} T") # A line of weight |c_n|^2 seen through a Gaussian window of length tau peaks at |c_n|^2 tau / sqrt(2 pi). ax3.vlines(res.energies, 0, fc * 0.8 * T_vib / np.sqrt(2 * np.pi), color="k", lw=1.5, label="exact Franck-Condon factors") center, width = V_morse(0.0, morse_params), abs(dV_morse(0.0, morse_params)) * sigma_x ax3.plot(E, np.exp(-0.5 * ((E - center) / width) ** 2) / (np.sqrt(2 * np.pi) * width), "k:", lw=1, label="Gaussian band from the slope") ax3.set_xlim(5, 35) ax3.set_xlabel("photon energy E") ax3.set_ylabel("absorption (wide lines: exact dynamics)") ax3.set_title("Bound upper surface: band and progression") ax3.legend(fontsize=7) for t_win, color in [(0.25 * T_vib, "0.5"), (0.8 * T_vib, "firebrick")]: ax4.plot(E, spectrum(E, times, C_rep_ex, t_win), color=color, lw=3, alpha=0.4) ax4.plot(E, spectrum(E, times, C_rep_tg, t_win), color=color, lw=1, ls="--", label=f"thawed, window {t_win / T_vib:.2f} T") center_r, width_r = V_rep(0.0, rep_params), abs(dV_rep(0.0, rep_params)) * sigma_x ax4.plot(E, np.exp(-0.5 * ((E - center_r) / width_r) ** 2) / (np.sqrt(2 * np.pi) * width_r), "k:", lw=1, label="Gaussian band from the slope") ax4.set_xlim(0, 35) ax4.set_xlabel("photon energy E") ax4.set_title("Repulsive upper surface: a structureless band") ax4.legend(fontsize=7) fig2.tight_layout() print(f"small-oscillation period on the upper surface T = {T_vib:.3f}, level spacing 2 pi hbar / T = {2 * np.pi / T_vib:.3f} at the bottom of the well") i_c = int(np.argmin(np.abs(res.energies - V_morse(0.0, morse_params)))) print(f"exact level spacing at the band centre: {res.energies[i_c + 1] - res.energies[i_c]:.3f}") for name, C_ex, C_tg in [("Morse", C_morse_ex, C_morse_tg), ("repulsive", C_rep_ex, C_rep_tg)]: for t_win in (0.25 * T_vib, 0.8 * T_vib): r = np.corrcoef(spectrum(E, times, C_ex, t_win), spectrum(E, times, C_tg, t_win))[0, 1] print(f"{name:9s}, window {t_win / T_vib:.2f} T: correlation of exact and thawed spectra {r:.3f}") plt.show() .. image-sg:: /api/gallery/semiclassical/propagators/images/sphx_glr_plot_heller_franck_condon_spectroscopy_002.png :alt: Bound upper surface: band and progression, Repulsive upper surface: a structureless band :srcset: /api/gallery/semiclassical/propagators/images/sphx_glr_plot_heller_franck_condon_spectroscopy_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none small-oscillation period on the upper surface T = 2.810, level spacing 2 pi hbar / T = 2.236 at the bottom of the well exact level spacing at the band centre: 1.674 Morse , window 0.25 T: correlation of exact and thawed spectra 0.999 Morse , window 0.80 T: correlation of exact and thawed spectra 0.850 repulsive, window 0.25 T: correlation of exact and thawed spectra 1.000 repulsive, window 0.80 T: correlation of exact and thawed spectra 1.000 .. GENERATED FROM PYTHON SOURCE LINES 212-217 Check ----- Short-time dynamics fix the band (the thawed Gaussian is exact there); resolving the Morse progression needs longer times, where the thawed Gaussian misses the anharmonicity; on the repulsive surface it stays exact. .. GENERATED FROM PYTHON SOURCE LINES 217-226 .. code-block:: Python corr = { (name, frac): np.corrcoef(spectrum(E, times, C_ex, frac * T_vib), spectrum(E, times, C_tg, frac * T_vib))[0, 1] for name, C_ex, C_tg in [("Morse", C_morse_ex, C_morse_tg), ("repulsive", C_rep_ex, C_rep_tg)] for frac in (0.25, 0.8) } assert corr["Morse", 0.25] > 0.99 and corr["repulsive", 0.25] > 0.99 and corr["repulsive", 0.8] > 0.99 assert corr["Morse", 0.8] < 0.95 assert abs(T_vib - 2 * np.pi / (a_m * np.sqrt(2 * D_e))) < 1e-12 assert res.energies[i_c + 1] - res.energies[i_c] < 2 * np.pi / T_vib .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 5.277 seconds) .. _sphx_glr_download_api_gallery_semiclassical_propagators_plot_heller_franck_condon_spectroscopy.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_heller_franck_condon_spectroscopy.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_heller_franck_condon_spectroscopy.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_heller_franck_condon_spectroscopy.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_