Computational Demo: Vibrational Transitions in CO Using Variational Method#
Jay Foley, University of North Carolina Charlotte
Summary#
In this activity, we study vibrational transitions of the carbon monoxide (CO) molecule using a quartic (fourth-order) polynomial approximation to the Morse potential. While this potential appears simpler than the Morse potential, the quartic oscillator energy eigenvalue equation cannot be solved analytically (even though the Morse oscillator can be). Nevertheless, it will give us some experience using approximate methods. In a related notebook, we computed perturbative corrections to vibrational energies using an anharmonic potential. Here, we will utilize the linear variational method to approximate the first fiew states of the quartic oscillator. Importantly, we will highlight the use of parity arguments to simplify computation of matrix elements in this approach.
# Quartic Oscillator Variational Method
# Using Harmonic Oscillator Basis
# Author: Jay Foley
# Purpose: Demonstrate how to construct the Hamiltonian matrix
# with analytic harmonic terms and numerical cubic + quartic perturbations (V')
import numpy as np
import matplotlib.pyplot as plt
from numpy import trapz
from scipy.special import hermite
from math import factorial
# ----------------------------------------------------
# 1. Constants and Physical Parameters (Atomic Units)
# ----------------------------------------------------
hbar_au = 1.0
k = 1.222147 # quadratic force constant
g = -4.462453 # cubic coefficient
h = 12.672998 # quartic coefficient
mu = 12506.108747 # reduced mass
r_eq = 2.132178 # equilibrium bond length
# ----------------------------------------------------
# 2. Helper functions (from user’s snippet, lightly adapted)
# ----------------------------------------------------
def compute_alpha(k, mu, hbar):
omega = np.sqrt(k / mu)
alpha = mu * omega / hbar
return alpha
def N(n, alpha):
return np.sqrt(1 / (2 ** n * factorial(n))) * (alpha / np.pi) ** 0.25
def psi(n, alpha, r, r_eq):
Hr = hermite(n)
psi_n = N(n, alpha) * Hr(np.sqrt(alpha) * (r - r_eq)) * np.exp(-0.5 * alpha * (r - r_eq)**2)
return psi_n
def harmonic_eigenvalue(n, k, mu, hbar):
"""
computes the harmonic oscillator eigenvalue for quantum number n
Arguments
---------
n : int
the quantum number
k : float
the force constant
mu : float
the reduced mass
hbar : float
reduced plancks constant
"""
return hbar * np.sqrt(k / mu) * (n + 0.5)
def potential_matrix_element(n, m, alpha, r, r_eq, V_p):
"""
computes the bra-ket <n|V_p|m> with basis functions \psi_n and \psi_m
Arguments
---------
n : int
the quantum number of the bra
m : int
the quantum number of the ket
alpha : float
mu * omega / hbar
r : np.array of floats
the position grid
r_eq : float
the equiibrium bond length
V_p : np.array of floats
the array of potential values along the r grid
"""
psi_n = psi(n, alpha, r, r_eq)
psi_m = psi(m, alpha, r, r_eq)
integrand = np.conj(psi_n) * V_p * psi_m
V_nm = np.trapz(integrand, r)
return V_nm
def build_hamiltonian_matrix(n_basis, k, alpha, r, r_eq, V_prime, mu, hbar):
H = np.zeros((n_basis, n_basis))
for n in range(n_basis):
for m in range(n_basis):
H0_nm = harmonic_eigenvalue(m, k, mu, hbar) * (n == m)
V_nm = potential_matrix_element(n, m, alpha, r, r_eq, V_prime)
H[n, m] = H0_nm + V_nm
return H
# ----------------------------------------------------
# 3. Construct potential and V'
# ----------------------------------------------------
# Grid around equilibrium
r_min = r_eq - 2.0
r_max = r_eq + 2.0
r_points = 2000
r = np.linspace(r_min, r_max, r_points)
# Full potential
V_total = 0.5 * k * (r - r_eq)**2 + (1/6)*g*(r - r_eq)**3 + (1/24)*h*(r - r_eq)**4
# Cubic part of potential V_cubic = (1/6)*g*x^3
V_cubic = (1/6)*g*(r - r_eq)**3
# Quartic part of the potential V_quartic = (1/24) * h * x^4
V_quartic = (1/24)*h*(r - r_eq)**4
# ----------------------------------------------------
# 4. Build and diagonalize Hamiltonian
# ----------------------------------------------------
alpha = compute_alpha(k, mu, hbar_au)
🧰 Helper Functions You’ll Need#
In this exercise, you’ll build each element of the \( 3 \times 3 \) Hamiltonian matrix by hand.
To do this, you will need to call two helper functions that are already defined for you.
⚙️ harmonic_eigenvalue(n, k, mu, hbar_au)#
Purpose:
Returns the energy eigenvalue of the unperturbed harmonic oscillator for a given quantum number \( n \):
Arguments:
n– integer, harmonic oscillator quantum numberk– harmonic force constant (in atomic units)mu– reduced mass (in atomic units)hbar_au– reduced Planck constant (1.0 in atomic units)
Returns:
A single float value for the harmonic energy $ E_n^{(0)} $ in atomic units.
Example:
# Compute the ground-state energy (n = 0)
E0 = harmonic_eigenvalue(0, k, mu, hbar_au)
print(E0)
⚙️ potential_matrix_element(n, m, alpha, r, r_eq, V_prime)#
Purpose:
Computes the matrix element
$\(
\langle n | V'(r) | m \rangle
\)\(
for a given potential \) V’(r) \(, using **numerical integration** over the position grid \) r $.
Here \(V'(r) \) can represent either:
the cubic perturbation:
$\( V_{\text{cubic}}(r) = \frac{1}{6} g (r - r_{\mathrm{eq}})^3 \)$or the quartic perturbation:
$\( V_{\text{quartic}}(r) = \frac{1}{24} h (r - r_{\mathrm{eq}})^4 \)$
Arguments:
n,m— integers
The quantum numbers for the harmonic oscillator bra and ket states.alpha— float
The width parameter \( \alpha = \frac{\mu \omega}{\hbar} \).r— NumPy array
The position grid (in atomic units) on which the integration is carried out.r_eq— float
The equilibrium bond length (in atomic units).V_prime— NumPy array
The perturbing potential \( V'(r) \), evaluated on the same grid.
Returns:
A single float value giving the matrix element
$\(
\langle n | V'(r) | m \rangle =
\int_{-\infty}^{\infty} \psi_n(r) \, V'(r) \, \psi_m(r) \, dr
\)$
The integration is performed numerically using the trapezoidal rule (np.trapz).
Example Usage:
# Example: compute ⟨0|V_cubic|1⟩
V_nm = potential_matrix_element(0, 1, alpha, r, r_eq, V_cubic)
print(V_nm)
# ----------------------------------------------------
# Guided Exercise: Build the 3x3 Hamiltonian Matrix by Hand
# ----------------------------------------------------
# Goal: Use the helper functions to fill in the Hamiltonian matrix elements
# for the cubic and quartic oscillator, in the harmonic oscillator basis.
# Reminder of selection rules:
# - Harmonic term: contributes ONLY to diagonal elements (analytic)
# - Cubic term: connects states of DIFFERENT parity (odd <-> even)
# - Quartic term: connects states of SAME parity (even <-> even, odd <-> odd)
# The potential terms:
# V_cubic = (1/6)*g*(r - r_eq)**3
# V_quartic = (1/24)*h*(r - r_eq)**4
# H[n, m] = <n|H0|m> + <n|V_cubic|m> + <n|V_quartic|m>
# ----------------------------------------------------
# STEP 1: Set matrix size and initialize
# ----------------------------------------------------
n_basis = 3
H_manual = np.zeros((n_basis, n_basis))
# ----------------------------------------------------
# STEP 2: Loop over n, m and fill in elements
# ----------------------------------------------------
for n in range(n_basis):
for m in range(n_basis):
# --- (a) Harmonic term: only diagonal ---
if n == m:
H0_nm = # insert call to harmonic_eigenvalue(n, k, mu, hbar_au)
else:
H0_nm = 0.0
# --- (b) Case 1: Bra and Ket have different parity ---
# Use potential_matrix_element() with V_cubic or V_quartic as appropriate
if (n + m) % 2 == 1: # This catches if n and m have different parity, so that the product of the bra and ket is an odd function!
Vo_nm = # insert call to potential_matrix_element(n, m, alpha, r, r_eq, ...)
else:
Vo_nm = 0.0
# --- (c) Case 2: Bra and Ket have same parity ---
# Use potential_matrix_element() with V_cubic or V_quartic as appropriate
if (n + m) % 2 == 0: # This catches if n and m have same parity, so that the product of the bra and ket is an even function!
Ve_nm = # insert call to potential_matrix_element(n, m, alpha, r, r_eq, ...)
else:
Ve_nm = 0.0
# --- (d) Combine all contributions ---
H_manual[n, m] = H0_nm + Vo_nm + Ve_nm
# ----------------------------------------------------
# STEP 3: Diagonalize and inspect results
# ----------------------------------------------------
E_manual, C_manual = np.linalg.eigh(H_manual)
# conversion from au to wn
au_to_invcm = 219474.63
print("3x3 Hamiltonian matrix (constructed manually):")
print(H_manual)
print("\nEigenvalues (in a.u.):")
for i, Ei in enumerate(E_manual):
print(f" State {i:2d}: E = {Ei:12.6f}")
print("\nEigenvalues (in cm^-1):")
for i, Ei in enumerate(E_manual):
print(f" State {i:2d}: E = {Ei * au_to_invcm:12.6f}")
# ----------------------------------------------------
# Display the variational coefficients for each state
# ----------------------------------------------------
# Each column of C_manual corresponds to an eigenvector (state)
# Each row corresponds to the coefficient for a particular basis state |nâź©
print("=== Variational Expansion Coefficients ===")
print("(Columns correspond to states; rows to basis functions)")
basis_labels = [f"|{n}âź©" for n in range(3)]
for state_index in range(3):
print(f"\nState {state_index}: Eigenvalue = {E_manual[state_index]* au_to_invcm:.6f} cm^-1")
print("Basis contributions:")
for n, label in enumerate(basis_labels):
coeff = C_manual[n, state_index]
coeff *= coeff
print(f" {label:<4} : {coeff: .6f}")