Vibrational Transitions in CO Using Perturbation Theory#
Jay Foley, University of North Carolina Charlotte
Learning Outcomes#
By the end of this workbook, students should be able to:
Explain the choice of zero-order Hamiltonian and perturbation for finding approximate solutions to the quartic oscillator Hamiltonian
Apply first- and second-order perturbation theory to obtain approximations to the fundamental transition
Compute perturbative corrections using the harmonic oscillator eigenstates and energies, and interpret how parity arguments can be used to simplify these computations
Convert between atomic unit and common unit systems
Assess the accuracy of perturbative approximations using the Morse oscillator as the ground truth
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 evaluating perturbative corrections to vibrational energies using an anharmonic potential. We can also compare exact solutions of the Morse and harmonic oscillator potentials with perturbative approximations. We will carry out the first-order energy corrections by hand, and utilize python to compute the second-order corrections. This provides an example of connecting quantum mechanical theory with computational practice.
Plot of model potentials#
# ——— Imports ———
import numpy as np
import matplotlib.pyplot as plt
# ——— Constants and conversions ———
De_eV = 11.225 # Dissociation energy (eV)
r_eq_ang = 1.1283 # Equilibrium bond length (Å)
mu_amu = 6.8606 # Reduced mass (amu)
beta_inv_ang = 2.30 # Morse curvature (Å⁻¹)
# Conversion factors
eV_to_au = 1 / 27.2114 # 1 Eh = 27.2114 eV
au_to_ang = 0.52917721067121 # 1 a0 = 0.529177 Å
amu_to_au = 1822.888486 # 1 amu = 1822.888 me
au_to_invcm = 219474.63 # Hartree to cm^-1
# Convert parameters to atomic units
De_au = De_eV * eV_to_au
r_eq_au = r_eq_ang / au_to_ang
mu_au = mu_amu * amu_to_au
beta_au = beta_inv_ang * au_to_ang
# ——— Helper function ———
def morse_potential(r, De, beta, r_eq):
"""
Evaluate the Morse potential at coordinate r (in atomic units).
Parameters
----------
r : array_like
Internuclear distance(s) in a0.
De : float
Dissociation energy (Hartrees).
beta : float
Range parameter (a0⁻¹).
r_eq : float
Equilibrium bond length (a0).
Returns
-------
V : ndarray
Morse potential in Hartrees.
"""
return De * (1 - np.exp(-beta * (r - r_eq)))**2
# ——— Define grid and potentials ———
r = np.linspace(0.5 * r_eq_au, 2.5 * r_eq_au, 500)
# Morse potential (exact)
V_Morse = morse_potential(r, De_au, beta_au, r_eq_au)
# Derivatives of Morse potential at r_eq (analytic)
k = 2 * De_au * beta_au**2
g = -6 * De_au * beta_au**3
h = 14 * De_au * beta_au**4
# Polynomial expansions around r_eq
x = r - r_eq_au
V_H = 0.5 * k * x**2 # pure harmonic
V_C = V_H + (1/6) * g * x**3 # harmonic + cubic
V_Q = V_C + (1/24) * h * x**4 # harmonic + cubic + quartic
# ——— Plotting ———
plt.style.use("seaborn-v0_8-colorblind")
fig, ax = plt.subplots(figsize=(7, 4.5))
ax.plot(r * au_to_ang, V_Morse * au_to_invcm, color="purple", lw=2.5, label="Morse potential ('ground truth')")
ax.plot(r * au_to_ang, V_H * au_to_invcm, "r--", lw=1.8, label="Harmonic (quadratic)")
ax.plot(r * au_to_ang, V_C * au_to_invcm, "b--", lw=1.8, label="Cubic expansion")
ax.plot(r * au_to_ang, V_Q * au_to_invcm, "g--", lw=1.8, label="Quartic expansion")
# Visual touches
ax.axvline(r_eq_au, color="gray", ls=":", lw=1)
ax.text(r_eq_au + 0.05, 0.02, r"$r_{eq}$", fontsize=11, color="gray")
ax.set_xlim(0.5 * r_eq_au * au_to_ang, 2.5 * r_eq_au * au_to_ang)
ax.set_ylim(0, De_au * 1.05 * au_to_invcm)
ax.set_xlabel("Bond length (Angstroms)", fontsize=12)
ax.set_ylabel("Potential Energy (cm$^{-1}$)", fontsize=12)
ax.set_title("Morse and Polynomial Approximations for CO", fontsize=13)
ax.legend(frameon=False, fontsize=10)
ax.grid(alpha=0.2)
plt.tight_layout()
plt.savefig("Quartic.png",dpi=300)
plt.show()
# print out the following constants in atomic units
print(f"Dissociation energy De: {De_au:.6f} a.u.")
print(f"Equilibrium bond length r_eq: {r_eq_au:.6f} a.u.")
print(f"Reduced mass mu: {mu_au:.6f} a.u.")
print(f"Morse parameter beta: {beta_au:.6f} a.u.")
# print out the force constants in atomic units
print(f"Force constant k: {k:.6f} a.u.")
print(f"Cubic constant g: {g:.6f} a.u.")
print(f"Quartic constant h: {h:.6f} a.u.")
Dissociation energy De: 0.412511 a.u.
Equilibrium bond length r_eq: 2.132178 a.u.
Reduced mass mu: 12506.108747 a.u.
Morse parameter beta: 1.217108 a.u.
Force constant k: 1.222147 a.u.
Cubic constant g: -4.462453 a.u.
Quartic constant h: 12.672998 a.u.
Full Quartic Oscillator Hamiltonian#
We can express the quartic oscillator Hamiltonian as $\( \hat{H} = -\frac{\hbar^2}{2\mu} \frac{d^2}{dx^2} + \frac{1}{2}k x^2 + \frac{1}{6} g x^3 + \frac{1}{24} h x^4. \)\( where here we understand \)x \equiv r - r_{eq}\( to be the bond length displaced by the equilibrium bond length such that the potential has a minumum at \)x = 0 \rightarrow r = r_{eq}$.
To draw some correspondence to the Morse oscillator model, we can write the Taylor expansion of the Morse potential around \(r-r_{eq}\)
\[
V_T(r) = \sum_{n=0}^\infty \frac{f^{(n)}(r_{eq})}{n!} (r - r_{eq})^n,
\]
where \(f^{(n)}(r_{eq})\) is the \(n^{th}\)-order derivative of the Morse potential evaluated at the equilibrium bondlength, e.g. \(f^{(1)}(r_{eq}) = \frac{d}{dr}V_{Morse}(r_{eq})\).
This leads to the following relations for each termin the quartic oscillator in terms of the Morse oscillator model
\[
k = \left. \frac{dV_{\text{Morse}}(r)}{dr} \right|_{r = r_{\text{eq}}} = 2 D_e \beta^2
\]
\[
g = \left. \frac{d^2V_{\text{Morse}}(r)}{dr^2} \right|_{r = r_{\text{eq}}}= -6 De \beta^3
\]
\[
h = \left. \frac{d^3V_{\text{Morse}}(r)}{dr^3} \right|_{r = r_{\text{eq}}} = 14 De \beta^4
\]
Harmonic Oscillator Reference#
We choose the harmonic oscillator as the zeroth-order Hamiltonian:
\[
\hat{H}_0 = -\frac{\hbar^2}{2\mu} \frac{d^2}{dx^2} + \frac{1}{2} k x^2, \qquad \omega=\sqrt{\frac{k}{\mu}},
\]
and define
\[
\alpha=\frac{\mu\omega}{\hbar}.
\]
The harmonic oscillator eigenfunctions are
\[
\psi_v^{(0)}(x)=\frac{1}{\sqrt{2^v v!}}\left(\frac{\alpha}{\pi}\right)^{1/4}
H_v(\sqrt{\alpha}\,x)\,e^{-\alpha x^2/2},
\]
where \(H_n\) are the physicists’ Hermite polynomials. A few relevant Hermite polynomials are
\(H_v(\sqrt{\alpha}\,x)\) |
parity |
|---|---|
\(H_0(\sqrt{\alpha}\,x) = 1\) |
even |
\(H_1(\sqrt{\alpha}\,x) = 2\sqrt{\alpha}\,x\) |
odd |
\(H_2(\sqrt{\alpha}\,x) = 4\alpha \,x^2-2\) |
even |
\(H_3(\sqrt{\alpha}\,x) = 8\alpha^{3/2} x^3 - 12 \sqrt{\alpha}\,x\) |
odd |
Question 1#
Why is
\[
E_v^{(0)} = \hbar \omega \left(v + \frac{1}{2}\right)
\]
an appropriate zeroth-order energy expression for the quartic oscillator?
a. This is an appropriate zeroth-order energy expression because \(v = 0 \) is the ground state
b. This is the expression for the Morse oscillator energy eigenvalues
c. This is the expression for the quartic oscillator energy eigenvalues
d. This is the expression for the exact energy eigenvalues of the unperturbed problem
Question 2#
Compute the zeroth-order approximation to the fundamental transition energy.
Express your answer in atomic units and in \(\,\text{cm}^{-1}\).
a. 0.004942 a.u., 1084 \(\text{cm}^{-1}\)
b. 0.01485 a.u., 3254 \(\text{cm}^{-1}\)
c. 0.004927 a.u., 1081 \(\text{cm}^{-1}\)
d. 0.009885 a.u., 2169 \(\text{cm}^{-1}\)
Question 3#
Why is
\[
V'(x) = \frac{g}{6}x^3 + \frac{h}{24}x^4
\]
an appropriate perturbation for the quartic oscillator model?
a. It contains the difference between the full potential and the potential of the zeroth-order Hamiltonian
b. This part of the potential makes the full Hamiltonian (harmonic + quartic) hard or impossible to solve analytically
c. This potential represents a relatively small correction to the zeroth-order potential (see plot of potentials)
d. All of the above
With this form of the perturbation, the first order corrections for the energy of a given state \(|\psi_v\rangle\) can be written as $\( E_v^{(1)} = \langle v|V'(x)|v\rangle = \int_{-\infty}^{\infty}\psi_v^{(0)}(x)\,V'(x)\,\psi_v^{(0)}(x)\,dx. \)\( Expanding the potential within these bra-kets gives: \)\( \frac{g}{6}\langle v|x^3|v\rangle + \frac{h}{24}\langle v|x^4|v\rangle = \frac{g}{6} \int_{-\infty}^{\infty}\psi_v^{(0)}(x)\,x^3\,\psi_v^{(0)}(x)\,dx + \frac{h}{24} \int_{-\infty}^{\infty}\psi_v^{(0)}(x)\,x^4\,\psi_v^{(0)}(x)\,dx. \)$
Question 4#
Using parity arguments, identify the bra-kets that must go to zero in the first-order corrections to any state
\(E_v^{(1)}\).
a. \(\langle v|x^4|v\rangle\) goes to zero for all \(v\)
b. \(\langle v|x^3|v\rangle\) goes to zero for all \(v\)
c. \(\langle v|x^4|v\rangle\) goes to zero for odd values of \(v\) (e.g. \(v = 1\)) and \(\langle v|x^3|v\rangle\) goes to zero for even values of \(v\) (e.g. \(v = 0\))
d. \(\langle v|x^4|v\rangle\) goes to zero for even values of \(v\) (e.g. \(v = 0\)) and \(\langle v|x^3|v\rangle\) goes to zero for odd values of \(v\) (e.g. \(v = 1\))
To address the last question, express the fundamental transition energy up to first order in perturbation theory. In particular, you need to compute the energy of the relevant states \(|\psi_{v_i}\rangle\) and \(|\psi_{v_f}\rangle\) that defines the fundamental transition up to first order and then compute the appropriate energy difference \(\Delta E_{fund} = E_{v_f} - E_{v_i}\) where each energy is approximated up to first order in perturbation theory.
Question 5 Which of the following expressions is the most accurate way to represent the fundamental transition energy for the quartic oscillator up to first order in perturbation theory?
a. \(\Delta E_{\text{fund}} = \big( E_1^{(0)} + E_1^{(1)} \big) - \big( E_0^{(0)} + E_0^{(1)} \big)\)
b. \(\Delta E_{\text{fund}} = \big( E_1^{(0)} - E_0^{(0)} \big) + E_1^{(1)}\)
c. \(\Delta E_{\text{fund}} = \big( E_1^{(0)} - E_0^{(0)} \big) + \big( E_1^{(1)} + E_0^{(1)} \big)\)
d. \(\Delta E_{\text{fund}} = E_1^{(0)} - E_0^{(0)}\)
Question 6 (for free response and discussion)#
Compute the first-order correction to the ground-state energy using
\[
E_0^{(1)} = \langle 0 | V'(x) | 0 \rangle = \int_{-\infty}^{\infty}\psi_v^{(0)}(x)\,V'(x)\,\psi_v^{(0)}(x)\,
\]
and express your answer in atomic units and in inverse centimeters.
Useful Integral Relation Definite integrals over all space of even-parity polynomials multiplied by Gaussians have the following form: $\( \int_{-\infty}^{\infty} x^{2n} {\rm e}^{-\alpha x^2} dx = \frac{\sqrt{\pi}}{2^n \alpha^{n + \frac{1}{2}}} \frac{(2n)!}{2^n n!} \)$
Question 7 (for free response and discussion)#
The first-order correction to the first excited-state energy is $\( E_1^{(1)} \approx 0.0001295 \; \text{a.u.} \approx 28.43 \: \text{cm}^{-1} \)$
Given this and your answer to Question 6, compute the fundamental transition energy of the quartic oscillator up to first order in perturbation theory and compare it to the fundamentals from the analogous harmonic oscillator and Morse oscillators:
Comparison of Fundamental Transition Energies#
Model |
Fundamental Frequency (cm⁻¹) |
|---|---|
Harmonic Oscillator |
2169. |
Morse Oscillator |
2143. |
Perturbation Theory (1st Order) |
— |
Second-order energy corrections#
We will now investigate how going to second order in perturbation theory will change our approximation to the fundamental transition. The code block below builds some helper functions to compute the requisite matrix elements, and then these are utilized in the subsequent code blocks to compute the fundamental up to second order.
from numpy import trapz
from scipy.special import hermite
from math import factorial
hbar_au = 1
def compute_alpha(k, mu, hbar):
""" Helper function to compute \alpha = \sqrt{k * \omega / \hbar}
Arguments
---------
k : float
the Harmonic force constant
mu : float
the reduced mass
hbar : float
reduced planck's constant
Returns
-------
alpha : float
\alpha = \sqrt{k * \omega / \hbar}
"""
# compute omega
omega = np.sqrt( k / mu )
# compute alpha
alpha = mu * omega / hbar
# return alpha
return alpha
def N(n, alpha):
""" Helper function to take the quantum number n of the Harmonic Oscillator and return the normalization constant
Arguments
---------
n : int
the quantum state of the harmonic oscillator
Returns
-------
N_n : float
the normalization constant
"""
return np.sqrt( 1 / (2 ** n * factorial(n)) ) * ( alpha / np.pi ) ** (1/4)
def psi(n, alpha, r, r_eq):
""" Helper function to evaluate the Harmonic Oscillator energy eigenfunction for state n
Arguments
---------
n : int
the quantum state of the harmonic oscillator
alpha : float
alpha value
r : float
position at which psi_n will be evaluated
r_eq : float
equilibrium bondlength
Returns
-------
psi_n : float
value of the harmonic oscillator energy eigenfunction
"""
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):
""" Helper function to evaluate the energy eigenvalue of the harmonic oscillator for state n"""
return hbar * np.sqrt(k/mu) * (n + 1/2)
def morse_eigenvalue(n, k, mu, De, hbar):
""" Helper function to evaluate the energy eigenvalue of the Morse oscillator for state n"""
omega = np.sqrt( k / mu )
xi = hbar * omega / (4 * De)
return hbar * omega * ( (n + 1/2) - xi * (n + 1/2) ** 2)
def potential_matrix_element(n, m, alpha, r, r_eq, V_p):
""" Helper function to compute <n|V_p|m> where V_p is the perturbing potential
Arguments
---------
n : int
quantum number of the bra state
m : int
quantum number of the ket state
alpha : float
alpha constant for bra/ket states
r : float
position grid for bra/ket states
r_eq : float
equilibrium bondlength for bra/ket states
V_p : float
potential array
Returns
-------
V_nm : float
<n | V_p | m >
"""
# bra
psi_n = psi(n, alpha, r, r_eq)
# ket
psi_m = psi(m, alpha, r, r_eq)
# integrand
integrand = np.conj(psi_n) * V_p * psi_m
# integrate
V_nm = np.trapz(integrand, r)
return V_nm
def build_hamiltonian_matrix(n_basis, k, alpha, r, r_eq, V_total, mu, hbar):
"""
Build the Hamiltonian matrix in the harmonic oscillator basis.
H = T + V, where matrix elements are:
H_nm = <n|T|m> + <n|V|m>
Arguments
---------
n_basis : int
number of basis functions to include (0, 1, ..., n_basis-1)
k : float
the harmonic foce constant
alpha : float
alpha constant for basis states
r : array
position grid
r_eq : float
equilibrium bond length
V_total : array
total potential energy along r grid
mu : float
reduced mass
hbar : float
reduced Planck's constant
Returns
-------
H : ndarray (n_basis x n_basis)
Hamiltonian matrix
"""
H = np.zeros((n_basis, n_basis))
for n in range(n_basis):
for m in range(n_basis):
# Kinetic energy contribution
T_nm = kinetic_matrix_element(n, m, alpha, mu, hbar)
# Potential energy contribution (using your existing function)
V_nm = potential_matrix_element(n, m, alpha, r, r_eq, V_total)
# Total Hamiltonian matrix element
H[n, m] = T_nm + V_nm
return H
def solve_variational(n_basis, alpha, r, r_eq, V_total, mu, hbar):
"""
Solve for eigenvalues and eigenvectors using the linear variational method.
Arguments
---------
n_basis : int
number of basis functions
alpha : float
alpha constant
r : array
position grid
r_eq : float
equilibrium bond length
V_total : array
total potential energy
mu : float
reduced mass
hbar : float
reduced Planck's constant
Returns
-------
eigenvalues : ndarray
energy eigenvalues (sorted)
eigenvectors : ndarray
corresponding eigenvectors (columns)
"""
# Build Hamiltonian matrix
H = build_hamiltonian_matrix(n_basis, alpha, r, r_eq, V_total, mu, hbar)
# Diagonalize
eigenvalues, eigenvectors = np.linalg.eigh(H)
return eigenvalues, eigenvectors
# compute alpha
alpha = compute_alpha(k, mu_au, hbar_au)
# compute psi_0 along the r grid
psi_0 =psi(0, alpha, r, r_eq_au)
# is it normalized?
Integral = trapz(psi_0 ** 2, r)
assert np.isclose(Integral, 1.0)
<>:6: SyntaxWarning: invalid escape sequence '\s'
<>:6: SyntaxWarning: invalid escape sequence '\s'
/var/folders/x2/ggtknn9s49s3dcq96ckcgfgh0000gn/T/ipykernel_43430/1733439826.py:6: SyntaxWarning: invalid escape sequence '\s'
""" Helper function to compute \alpha = \sqrt{k * \omega / \hbar}
/var/folders/x2/ggtknn9s49s3dcq96ckcgfgh0000gn/T/ipykernel_43430/1733439826.py:144: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.
Integral = trapz(psi_0 ** 2, r)
E0_harmonic = harmonic_eigenvalue(0, k, mu_au, hbar_au)
E0_morse = morse_eigenvalue(0, k, mu_au, De_au, hbar_au)
E1_harmonic = harmonic_eigenvalue(1, k, mu_au, hbar_au)
E1_morse = morse_eigenvalue(1, k, mu_au, De_au, hbar_au)
# define the perturbing potential as the the difference between V_Q and V_H
V_pert = V_Q - V_H
# compute first order corrections
E0_fo_correction = potential_matrix_element(0, 0, alpha, r, r_eq_au, V_pert)
d = h / 24
_expected_E0_fo = ( 3 * d) / (4 * alpha ** 2)
E1_fo_correction = potential_matrix_element(1, 1, alpha, r, r_eq_au, V_pert)
# compute energy up to first order
E1_pt1 = E1_harmonic + E1_fo_correction
E0_pt1 = E0_harmonic + E0_fo_correction
# compute fundamental frequencies at different levels of theory
fundamental_harmonic_au = E1_harmonic - E0_harmonic
fundamental_morse_au = E1_morse - E0_morse
fundamental_pt1_au = E1_pt1 - E0_pt1
# print results
print(F"First order correction to ground state energy: {E0_fo_correction * au_to_invcm} cm^-1")
print(F"Expected first order correction to ground state energy: {_expected_E0_fo * au_to_invcm} cm^-1")
print(F"Difference between computed and expected: {(E0_fo_correction - _expected_E0_fo) * au_to_invcm} cm^-1")
print(F"First order correction to first excited state energy: {E1_fo_correction} au")
print(F"First order correction to first excited state energy: {E1_fo_correction * au_to_invcm} cm^-1")
print("Fundamental Frequency (Harmonic) : ", fundamental_harmonic_au * au_to_invcm, " cm^-1")
print("Fundeamental Frequency in au (Harmonic): ", fundamental_harmonic_au)
print("Fundamental Frequency (Morse) : ", fundamental_morse_au * au_to_invcm, " cm^-1")
print("Fundamental Frequency (PT1) : ", fundamental_pt1_au * au_to_invcm, " cm^-1")
First order correction to ground state energy: 5.686802216774283 cm^-1
Expected first order correction to ground state energy: 5.686802216774276 cm^-1
Difference between computed and expected: 6.692480737072095e-15 cm^-1
First order correction to first excited state energy: 0.00012955488788782298 au
First order correction to first excited state energy: 28.434011083871432 cm^-1
Fundamental Frequency (Harmonic) : 2169.626233701669 cm^-1
Fundeamental Frequency in au (Harmonic): 0.009885544555658526
Fundamental Frequency (Morse) : 2143.6294235678447 cm^-1
Fundamental Frequency (PT1) : 2192.3734425687667 cm^-1
/var/folders/x2/ggtknn9s49s3dcq96ckcgfgh0000gn/T/ipykernel_43430/1733439826.py:133: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.
V_nm = np.trapz(integrand, r)
# set up 2nd order correction to ground state energy and first excited state energy
E0_so_correction = 0.0
E1_so_correction = 0.0
for m in range(0, 10):
if m != 0:
V_0m = potential_matrix_element(0, m, alpha, r, r_eq_au, V_pert)
E0_so_correction += ( np.abs(V_0m) ** 2 ) / ( E0_harmonic - harmonic_eigenvalue(m, k, mu_au, hbar_au) )
print(F"m={m}, V_0m = {V_0m * au_to_invcm} cm^-1, contribution to E0_so: {( np.abs(V_0m) ** 2 ) / ( E0_harmonic - harmonic_eigenvalue(m, k, mu_au, hbar_au) ) * au_to_invcm} cm^-1")
if m != 1:
V_1m = potential_matrix_element(1, m, alpha, r, r_eq_au, V_pert)
E1_so_correction += ( np.abs(V_1m) ** 2 ) / ( E1_harmonic - harmonic_eigenvalue(m, k, mu_au, hbar_au) )
print(F"m={m}, V_1m = {V_1m * au_to_invcm} cm^-1, contribution to E1_so: {( np.abs(V_1m) ** 2 ) / ( E1_harmonic - harmonic_eigenvalue(m, k, mu_au, hbar_au) ) * au_to_invcm} cm^-1")
E0_pt2 = E0_pt1 + E0_so_correction
E1_pt2 = E1_pt1 + E1_so_correction
fundamental_pt2_au = E1_pt2 - E0_pt2
print(F"Second order correction to ground state energy: {E0_so_correction * au_to_invcm} cm^-1")
print(F"Second order correction to first excited state energy: {E1_so_correction * au_to_invcm} cm^-1")
print("Fundamental Frequency (PT2) : ", fundamental_pt2_au * au_to_invcm, " cm^-1")
m=0, V_1m = -125.9501701232201 cm^-1, contribution to E1_so: 7.311602850138361 cm^-1
m=1, V_0m = -125.9501701232201 cm^-1, contribution to E0_so: -7.311602850138361 cm^-1
m=2, V_0m = 16.084705642991135 cm^-1, contribution to E0_so: -0.0596226556452224 cm^-1
m=2, V_1m = -356.2408775429131 cm^-1, contribution to E1_so: -58.492822801106925 cm^-1
m=3, V_0m = -102.83788327287469 cm^-1, contribution to E0_so: -1.6248006333640819 cm^-1
m=3, V_1m = 46.43254566408416 cm^-1, contribution to E1_so: -0.4968554637101873 cm^-1
m=4, V_0m = 9.286509132816835 cm^-1, contribution to E0_so: -0.009937109274203753 cm^-1
m=4, V_1m = -205.67576654574958 cm^-1, contribution to E1_so: -6.499202533456339 cm^-1
m=5, V_0m = -6.023232663364886e-14 cm^-1, contribution to E0_so: -3.3442932384835987e-31 cm^-1
m=5, V_1m = 20.76526569465107 cm^-1, contribution to E1_so: -0.04968554637101877 cm^-1
m=6, V_0m = -7.436089707857883e-15 cm^-1, contribution to E0_so: -4.2476924738144574e-33 cm^-1
m=6, V_1m = -1.1228495458865403e-13 cm^-1, contribution to E1_so: -1.1622196331453212e-30 cm^-1
m=7, V_0m = 1.4988368317401045e-14 cm^-1, contribution to E0_so: -1.4791960893571368e-32 cm^-1
m=7, V_1m = -3.494962162693205e-14 cm^-1, contribution to E1_so: -9.383152674656135e-32 cm^-1
m=8, V_0m = -8.737405406733012e-15 cm^-1, contribution to E0_so: -4.398352816245065e-33 cm^-1
m=8, V_1m = 2.2494171366270097e-14 cm^-1, contribution to E1_so: -3.331632080845579e-32 cm^-1
m=9, V_0m = -5.3446894775228535e-15 cm^-1, contribution to E0_so: -1.4629097126613594e-33 cm^-1
m=9, V_1m = 1.5615788386501556e-14 cm^-1, contribution to E1_so: -1.404924285714132e-32 cm^-1
Second order correction to ground state energy: -9.00596324842187 cm^-1
Second order correction to first excited state energy: -58.22696349450612 cm^-1
Fundamental Frequency (PT2) : 2143.1524423226824 cm^-1
/var/folders/x2/ggtknn9s49s3dcq96ckcgfgh0000gn/T/ipykernel_43430/1733439826.py:133: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.
V_nm = np.trapz(integrand, r)
# ============================================================================
# Apply variational method to quartic oscillator
# ============================================================================
print("\n" + "="*70)
print("LINEAR VARIATIONAL METHOD RESULTS")
print("="*70)
# Test with different basis set sizes
for n_basis in [5, 10, 15, 20]:
# Solve using variational method with V_Q (quartic potential)
eigenvalues, eigenvectors = solve_variational(
n_basis, alpha, r, r_eq_au, V_Q, mu_au, hbar_au
)
E0_var = eigenvalues[0]
E1_var = eigenvalues[1]
fundamental_var_au = E1_var - E0_var
print(f"\nBasis size: {n_basis}")
print(f" E0 (variational) : {E0_var * au_to_invcm:.4f} cm^-1")
print(f" E1 (variational) : {E1_var * au_to_invcm:.4f} cm^-1")
print(f" Fundamental (var) : {fundamental_var_au * au_to_invcm:.4f} cm^-1")
# Use largest basis for comparison
n_basis = 20
eigenvalues, eigenvectors = solve_variational(
n_basis, alpha, r, r_eq_au, V_Q, mu_au, hbar_au
)
E0_var = eigenvalues[0]
E1_var = eigenvalues[1]
fundamental_var_au = E1_var - E0_var
print("\n" + "="*70)
print("COMPARISON OF METHODS (n_basis = 20)")
print("="*70)
print(f"Ground state energy:")
print(f" Harmonic : {E0_harmonic * au_to_invcm:.4f} cm^-1")
print(f" PT1 : {E0_pt1 * au_to_invcm:.4f} cm^-1")
print(f" PT2 : {E0_pt2 * au_to_invcm:.4f} cm^-1")
print(f" Variational : {E0_var * au_to_invcm:.4f} cm^-1")
print(f" Morse (exact) : {E0_morse * au_to_invcm:.4f} cm^-1")
print(f"\nFirst excited state energy:")
print(f" Harmonic : {E1_harmonic * au_to_invcm:.4f} cm^-1")
print(f" PT1 : {E1_pt1 * au_to_invcm:.4f} cm^-1")
print(f" PT2 : {E1_pt2 * au_to_invcm:.4f} cm^-1")
print(f" Variational : {E1_var * au_to_invcm:.4f} cm^-1")
print(f" Morse (exact) : {E1_morse * au_to_invcm:.4f} cm^-1")
print(f"\nFundamental transition:")
print(f" Harmonic : {fundamental_harmonic_au * au_to_invcm:.4f} cm^-1")
print(f" PT1 : {fundamental_pt1_au * au_to_invcm:.4f} cm^-1")
print(f" PT2 : {fundamental_pt2_au * au_to_invcm:.4f} cm^-1")
print(f" Variational : {fundamental_var_au * au_to_invcm:.4f} cm^-1")
print(f" Morse (exact) : {fundamental_morse_au * au_to_invcm:.4f} cm^-1")
print(f"\nError in fundamental vs Morse:")
print(f" Harmonic : {(fundamental_harmonic_au - fundamental_morse_au) * au_to_invcm:.4f} cm^-1")
print(f" PT1 : {(fundamental_pt1_au - fundamental_morse_au) * au_to_invcm:.4f} cm^-1")
print(f" PT2 : {(fundamental_pt2_au - fundamental_morse_au) * au_to_invcm:.4f} cm^-1")
print(f" Variational : {(fundamental_var_au - fundamental_morse_au) * au_to_invcm:.4f} cm^-1")
# Optional: Print expansion coefficients for ground state
print(f"\nGround state expansion coefficients (first 10):")
for i in range(min(10, n_basis)):
print(f" c_{i} = {eigenvectors[i, 0]:.6f}")
======================================================================
LINEAR VARIATIONAL METHOD RESULTS
======================================================================
Basis size: 5
E0 (variational) : 1529.5581 cm^-1
E1 (variational) : 4588.3054 cm^-1
Fundamental (var) : 3058.7473 cm^-1
Basis size: 10
E0 (variational) : 1528.0493 cm^-1
E1 (variational) : 4548.2488 cm^-1
Fundamental (var) : 3020.1995 cm^-1
Basis size: 15
E0 (variational) : 1528.0451 cm^-1
E1 (variational) : 4547.9982 cm^-1
Fundamental (var) : 3019.9530 cm^-1
Basis size: 20
E0 (variational) : 1528.0451 cm^-1
E1 (variational) : 4547.9969 cm^-1
Fundamental (var) : 3019.9518 cm^-1
======================================================================
COMPARISON OF METHODS (n_basis = 20)
======================================================================
Ground state energy:
Harmonic : 1084.8131 cm^-1
PT1 : 1090.4999 cm^-1
PT2 : 1081.4940 cm^-1
Variational : 1528.0451 cm^-1
Morse (exact) : 1081.5635 cm^-1
First excited state energy:
Harmonic : 3254.4394 cm^-1
PT1 : 3282.8734 cm^-1
PT2 : 3224.6464 cm^-1
Variational : 4547.9969 cm^-1
Morse (exact) : 3225.1929 cm^-1
Fundamental transition:
Harmonic : 2169.6262 cm^-1
PT1 : 2192.3734 cm^-1
PT2 : 2143.1524 cm^-1
Variational : 3019.9518 cm^-1
Morse (exact) : 2143.6294 cm^-1
Error in fundamental vs Morse:
Harmonic : 25.9968 cm^-1
PT1 : 48.7440 cm^-1
PT2 : -0.4770 cm^-1
Variational : 876.3224 cm^-1
Ground state expansion coefficients (first 10):
c_0 = 0.989566
c_1 = 0.062785
c_2 = 0.123691
c_3 = 0.030783
c_4 = 0.021239
c_5 = 0.009367
c_6 = 0.004791
c_7 = 0.002511
c_8 = 0.001269
c_9 = 0.000674
/var/folders/x2/ggtknn9s49s3dcq96ckcgfgh0000gn/T/ipykernel_43430/1733439826.py:133: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.
V_nm = np.trapz(integrand, r)
# ============================================================================
# PERTURBATION THEORY CALCULATIONS
# ============================================================================
# Compute reference energies
E0_harmonic = harmonic_eigenvalue(0, k, mu_au, hbar_au)
E1_harmonic = harmonic_eigenvalue(1, k, mu_au, hbar_au)
E0_morse = morse_eigenvalue(0, k, mu_au, De_au, hbar_au)
E1_morse = morse_eigenvalue(1, k, mu_au, De_au, hbar_au)
# Define the perturbing potential as the difference between V_Q and V_H
V_pert = V_Q - V_H
print("\n" + "="*70)
print("FIRST ORDER PERTURBATION THEORY")
print("="*70)
# Compute first order corrections
E0_fo_correction = potential_matrix_element(0, 0, alpha, r, r_eq_au, V_pert)
E1_fo_correction = potential_matrix_element(1, 1, alpha, r, r_eq_au, V_pert)
# Analytical expectation for ground state (for validation)
d = h / 24
_expected_E0_fo = (3 * d) / (4 * alpha ** 2)
print(f"\nGround state (n=0):")
print(f" First order correction : {E0_fo_correction * au_to_invcm:>10.4f} cm^-1")
print(f" Expected (analytical) : {_expected_E0_fo * au_to_invcm:>10.4f} cm^-1")
print(f" Difference : {(E0_fo_correction - _expected_E0_fo) * au_to_invcm:>10.6f} cm^-1")
print(f"\nFirst excited state (n=1):")
print(f" First order correction : {E1_fo_correction * au_to_invcm:>10.4f} cm^-1")
# Compute energies up to first order
E0_pt1 = E0_harmonic + E0_fo_correction
E1_pt1 = E1_harmonic + E1_fo_correction
fundamental_pt1_au = E1_pt1 - E0_pt1
print(f"\nTotal energies (PT1):")
print(f" E0 (PT1) : {E0_pt1 * au_to_invcm:>10.4f} cm^-1")
print(f" E1 (PT1) : {E1_pt1 * au_to_invcm:>10.4f} cm^-1")
print(f" Fundamental transition (PT1) : {fundamental_pt1_au * au_to_invcm:>10.4f} cm^-1")
# ============================================================================
print("\n" + "="*70)
print("SECOND ORDER PERTURBATION THEORY")
print("="*70)
# Set up 2nd order correction to ground state and first excited state energy
E0_so_correction = 0.0
E1_so_correction = 0.0
n_states = 10 # number of states to sum over
print(f"\nSumming over {n_states} intermediate states...\n")
# Detailed contributions for ground state
print("Ground state (n=0) contributions:")
for m in range(n_states):
if m != 0:
V_0m = potential_matrix_element(0, m, alpha, r, r_eq_au, V_pert)
Em = harmonic_eigenvalue(m, k, mu_au, hbar_au)
contribution = (np.abs(V_0m) ** 2) / (E0_harmonic - Em)
E0_so_correction += contribution
print(f" m={m:2d} |V_0m|² = {np.abs(V_0m)**2 * au_to_invcm**2:>12.4e} (cm^-1)² "
f"→ ΔE = {contribution * au_to_invcm:>10.6f} cm^-1")
print(f"\nFirst excited state (n=1) contributions:")
for m in range(n_states):
if m != 1:
V_1m = potential_matrix_element(1, m, alpha, r, r_eq_au, V_pert)
Em = harmonic_eigenvalue(m, k, mu_au, hbar_au)
contribution = (np.abs(V_1m) ** 2) / (E1_harmonic - Em)
E1_so_correction += contribution
print(f" m={m:2d} |V_1m|² = {np.abs(V_1m)**2 * au_to_invcm**2:>12.4e} (cm^-1)² "
f"→ ΔE = {contribution * au_to_invcm:>10.6f} cm^-1")
print(f"\nSecond order corrections:")
print(f" Ground state (n=0) : {E0_so_correction * au_to_invcm:>10.4f} cm^-1")
print(f" First excited state (n=1) : {E1_so_correction * au_to_invcm:>10.4f} cm^-1")
# Compute energies up to second order
E0_pt2 = E0_pt1 + E0_so_correction
E1_pt2 = E1_pt1 + E1_so_correction
fundamental_pt2_au = E1_pt2 - E0_pt2
print(f"\nTotal energies (PT2):")
print(f" E0 (PT2) : {E0_pt2 * au_to_invcm:>10.4f} cm^-1")
print(f" E1 (PT2) : {E1_pt2 * au_to_invcm:>10.4f} cm^-1")
print(f" Fundamental transition (PT2) : {fundamental_pt2_au * au_to_invcm:>10.4f} cm^-1")
# ============================================================================
print("\n" + "="*70)
print("SUMMARY: COMPARISON OF ALL METHODS")
print("="*70)
fundamental_harmonic_au = E1_harmonic - E0_harmonic
fundamental_morse_au = E1_morse - E0_morse
print(f"\nFundamental transition frequency:")
print(f" Harmonic : {fundamental_harmonic_au * au_to_invcm:>10.4f} cm^-1")
print(f" PT1 : {fundamental_pt1_au * au_to_invcm:>10.4f} cm^-1")
print(f" PT2 : {fundamental_pt2_au * au_to_invcm:>10.4f} cm^-1")
print(f" Morse (reference) : {fundamental_morse_au * au_to_invcm:>10.4f} cm^-1")
print(f"\nError vs Morse:")
print(f" Harmonic : {(fundamental_harmonic_au - fundamental_morse_au) * au_to_invcm:>+10.4f} cm^-1")
print(f" PT1 : {(fundamental_pt1_au - fundamental_morse_au) * au_to_invcm:>+10.4f} cm^-1")
print(f" PT2 : {(fundamental_pt2_au - fundamental_morse_au) * au_to_invcm:>+10.4f} cm^-1")
print("\n" + "="*70)
======================================================================
FIRST ORDER PERTURBATION THEORY
======================================================================
Ground state (n=0):
First order correction : 5.6868 cm^-1
Expected (analytical) : 5.6868 cm^-1
Difference : 0.000000 cm^-1
First excited state (n=1):
First order correction : 28.4340 cm^-1
Total energies (PT1):
E0 (PT1) : 1090.4999 cm^-1
E1 (PT1) : 3282.8734 cm^-1
Fundamental transition (PT1) : 2192.3734 cm^-1
======================================================================
SECOND ORDER PERTURBATION THEORY
======================================================================
Summing over 10 intermediate states...
Ground state (n=0) contributions:
m= 1 |V_0m|² = 1.5863e+04 (cm^-1)² → ΔE = -7.311603 cm^-1
m= 2 |V_0m|² = 2.5872e+02 (cm^-1)² → ΔE = -0.059623 cm^-1
m= 3 |V_0m|² = 1.0576e+04 (cm^-1)² → ΔE = -1.624801 cm^-1
m= 4 |V_0m|² = 8.6239e+01 (cm^-1)² → ΔE = -0.009937 cm^-1
m= 5 |V_0m|² = 3.6279e-27 (cm^-1)² → ΔE = -0.000000 cm^-1
m= 6 |V_0m|² = 5.5295e-29 (cm^-1)² → ΔE = -0.000000 cm^-1
m= 7 |V_0m|² = 2.2465e-28 (cm^-1)² → ΔE = -0.000000 cm^-1
m= 8 |V_0m|² = 7.6342e-29 (cm^-1)² → ΔE = -0.000000 cm^-1
m= 9 |V_0m|² = 2.8566e-29 (cm^-1)² → ΔE = -0.000000 cm^-1
First excited state (n=1) contributions:
m= 0 |V_1m|² = 1.5863e+04 (cm^-1)² → ΔE = 7.311603 cm^-1
m= 2 |V_1m|² = 1.2691e+05 (cm^-1)² → ΔE = -58.492823 cm^-1
m= 3 |V_1m|² = 2.1560e+03 (cm^-1)² → ΔE = -0.496855 cm^-1
m= 4 |V_1m|² = 4.2303e+04 (cm^-1)² → ΔE = -6.499203 cm^-1
m= 5 |V_1m|² = 4.3120e+02 (cm^-1)² → ΔE = -0.049686 cm^-1
m= 6 |V_1m|² = 1.2608e-26 (cm^-1)² → ΔE = -0.000000 cm^-1
m= 7 |V_1m|² = 1.2215e-27 (cm^-1)² → ΔE = -0.000000 cm^-1
m= 8 |V_1m|² = 5.0599e-28 (cm^-1)² → ΔE = -0.000000 cm^-1
m= 9 |V_1m|² = 2.4385e-28 (cm^-1)² → ΔE = -0.000000 cm^-1
Second order corrections:
Ground state (n=0) : -9.0060 cm^-1
First excited state (n=1) : -58.2270 cm^-1
Total energies (PT2):
E0 (PT2) : 1081.4940 cm^-1
E1 (PT2) : 3224.6464 cm^-1
Fundamental transition (PT2) : 2143.1524 cm^-1
======================================================================
SUMMARY: COMPARISON OF ALL METHODS
======================================================================
Fundamental transition frequency:
Harmonic : 2169.6262 cm^-1
PT1 : 2192.3734 cm^-1
PT2 : 2143.1524 cm^-1
Morse (reference) : 2143.6294 cm^-1
Error vs Morse:
Harmonic : +25.9968 cm^-1
PT1 : +48.7440 cm^-1
PT2 : -0.4770 cm^-1
======================================================================
/var/folders/x2/ggtknn9s49s3dcq96ckcgfgh0000gn/T/ipykernel_43430/1733439826.py:133: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.
V_nm = np.trapz(integrand, r)
def build_kinetic_matrix_fd(r, mu, hbar):
"""
Build kinetic energy matrix using 5-point stencil finite difference.
The kinetic energy operator is:
T = -hbar^2 / (2*mu) * d^2/dr^2
5-point stencil for second derivative:
d^2f/dr^2|_i ≈ (-f_{i-2} + 16*f_{i-1} - 30*f_i + 16*f_{i+1} - f_{i+2}) / (12*h^2)
Arguments
---------
r : array
position grid (equally spaced)
mu : float
reduced mass
hbar : float
reduced Planck's constant
Returns
-------
T : ndarray (n_grid x n_grid)
Kinetic energy matrix
"""
n_grid = len(r)
h = r[1] - r[0] # grid spacing (assumes uniform grid)
# Initialize kinetic energy matrix
T = np.zeros((n_grid, n_grid))
# Prefactor for kinetic energy
prefactor = -hbar**2 / (2 * mu * 12 * h**2)
# Fill in the matrix using 5-point stencil
# Note: We use boundary conditions where wavefunction goes to zero at edges
for i in range(n_grid):
# Diagonal term (coefficient of psi_i)
T[i, i] = prefactor * (-30)
# Off-diagonal terms
if i >= 1:
T[i, i-1] = prefactor * 16 # coefficient of psi_{i-1}
if i <= n_grid - 2:
T[i, i+1] = prefactor * 16 # coefficient of psi_{i+1}
if i >= 2:
T[i, i-2] = prefactor * (-1) # coefficient of psi_{i-2}
if i <= n_grid - 3:
T[i, i+2] = prefactor * (-1) # coefficient of psi_{i+2}
return T
def build_potential_matrix_fd(V):
"""
Build potential energy matrix (diagonal).
Arguments
---------
V : array
potential energy values along the grid
Returns
-------
V_mat : ndarray (n_grid x n_grid)
Potential energy matrix (diagonal)
"""
return np.diag(V)
def solve_finite_difference(r, V, mu, hbar, n_states=10):
"""
Solve Schrödinger equation using finite difference method.
Arguments
---------
r : array
position grid
V : array
potential energy along grid
mu : float
reduced mass
hbar : float
reduced Planck's constant
n_states : int
number of lowest eigenvalues/eigenvectors to return
Returns
-------
eigenvalues : ndarray
energy eigenvalues (sorted, lowest n_states)
eigenvectors : ndarray
corresponding eigenvectors (columns, normalized)
"""
# Build kinetic and potential matrices
T = build_kinetic_matrix_fd(r, mu, hbar)
V_mat = build_potential_matrix_fd(V)
# Total Hamiltonian
H = T + V_mat
# Diagonalize (eigh for symmetric/Hermitian matrices)
eigenvalues, eigenvectors = np.linalg.eigh(H)
# Normalize eigenvectors properly (account for grid spacing)
h = r[1] - r[0]
for i in range(eigenvectors.shape[1]):
norm = np.sqrt(np.trapz(eigenvectors[:, i]**2, r))
eigenvectors[:, i] /= norm
# Return only the lowest n_states
return eigenvalues[:n_states], eigenvectors[:, :n_states]
# ============================================================================
# Apply finite difference method to quartic oscillator
# ============================================================================
print("\n" + "="*70)
print("FINITE DIFFERENCE METHOD RESULTS")
print("="*70)
# Solve using finite difference with V_Q (quartic potential)
eigenvalues_fd, eigenvectors_fd = solve_finite_difference(
r, V_Q, mu_au, hbar_au, n_states=10
)
E0_fd = eigenvalues_fd[0]
E1_fd = eigenvalues_fd[1]
fundamental_fd_au = E1_fd - E0_fd
print(f"\nFinite Difference Results:")
print(f" E0 (FD) : {E0_fd * au_to_invcm:>10.4f} cm^-1")
print(f" E1 (FD) : {E1_fd * au_to_invcm:>10.4f} cm^-1")
print(f" Fundamental (FD) : {fundamental_fd_au * au_to_invcm:>10.4f} cm^-1")
print(f"\nFirst 5 eigenvalues (cm^-1):")
for i in range(5):
print(f" E{i} = {eigenvalues_fd[i] * au_to_invcm:>10.4f} cm^-1")
# ============================================================================
# Compare all methods
# ============================================================================
print("\n" + "="*70)
print("COMPREHENSIVE COMPARISON OF ALL METHODS")
print("="*70)
# Variational with N=20 for comparison
eigenvalues_var, eigenvectors_var = solve_variational(
20, alpha, r, r_eq_au, V_Q, mu_au, hbar_au
)
E0_var = eigenvalues_var[0]
E1_var = eigenvalues_var[1]
fundamental_var_au = E1_var - E0_var
print(f"\nGround state energy (E0):")
print(f" Harmonic : {E0_harmonic * au_to_invcm:>10.4f} cm^-1")
print(f" PT1 : {E0_pt1 * au_to_invcm:>10.4f} cm^-1")
print(f" PT2 : {E0_pt2 * au_to_invcm:>10.4f} cm^-1")
print(f" Variational (N=20) : {E0_var * au_to_invcm:>10.4f} cm^-1")
print(f" Finite Difference : {E0_fd * au_to_invcm:>10.4f} cm^-1")
print(f" Morse (reference) : {E0_morse * au_to_invcm:>10.4f} cm^-1")
print(f"\nFirst excited state energy (E1):")
print(f" Harmonic : {E1_harmonic * au_to_invcm:>10.4f} cm^-1")
print(f" PT1 : {E1_pt1 * au_to_invcm:>10.4f} cm^-1")
print(f" PT2 : {E1_pt2 * au_to_invcm:>10.4f} cm^-1")
print(f" Variational (N=20) : {E1_var * au_to_invcm:>10.4f} cm^-1")
print(f" Finite Difference : {E1_fd * au_to_invcm:>10.4f} cm^-1")
print(f" Morse (reference) : {E1_morse * au_to_invcm:>10.4f} cm^-1")
print(f"\nFundamental transition:")
print(f" Harmonic : {fundamental_harmonic_au * au_to_invcm:>10.4f} cm^-1")
print(f" PT1 : {fundamental_pt1_au * au_to_invcm:>10.4f} cm^-1")
print(f" PT2 : {fundamental_pt2_au * au_to_invcm:>10.4f} cm^-1")
print(f" Variational (N=20) : {fundamental_var_au * au_to_invcm:>10.4f} cm^-1")
print(f" Finite Difference : {fundamental_fd_au * au_to_invcm:>10.4f} cm^-1")
print(f" Morse (reference) : {fundamental_morse_au * au_to_invcm:>10.4f} cm^-1")
print(f"\nError in fundamental vs Morse:")
print(f" Harmonic : {(fundamental_harmonic_au - fundamental_morse_au) * au_to_invcm:>+10.4f} cm^-1")
print(f" PT1 : {(fundamental_pt1_au - fundamental_morse_au) * au_to_invcm:>+10.4f} cm^-1")
print(f" PT2 : {(fundamental_pt2_au - fundamental_morse_au) * au_to_invcm:>+10.4f} cm^-1")
print(f" Variational (N=20) : {(fundamental_var_au - fundamental_morse_au) * au_to_invcm:>+10.4f} cm^-1")
print(f" Finite Difference : {(fundamental_fd_au - fundamental_morse_au) * au_to_invcm:>+10.4f} cm^-1")
print("\n" + "="*70)
# Optional: Plot wavefunctions
import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
# Plot ground state wavefunction
ax1.plot(r * au_to_ang, eigenvectors_fd[:, 0], 'b-', linewidth=2, label='FD: n=0')
ax1.plot(r * au_to_ang, psi(0, alpha, r, r_eq_au), 'r--', linewidth=2, label='Harmonic: n=0')
ax1.set_xlabel('r (Å)', fontsize=12)
ax1.set_ylabel('ψ₀(r)', fontsize=12)
ax1.set_title('Ground State Wavefunction', fontsize=14)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
# Plot first excited state wavefunction
ax2.plot(r * au_to_ang, eigenvectors_fd[:, 1], 'b-', linewidth=2, label='FD: n=1')
ax2.plot(r * au_to_ang, psi(1, alpha, r, r_eq_au), 'r--', linewidth=2, label='Harmonic: n=1')
ax2.set_xlabel('r (Å)', fontsize=12)
ax2.set_ylabel('ψ₁(r)', fontsize=12)
ax2.set_title('First Excited State Wavefunction', fontsize=14)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print("\nWavefunction comparison plots generated!")
======================================================================
FINITE DIFFERENCE METHOD RESULTS
======================================================================
Finite Difference Results:
E0 (FD) : 1081.6974 cm^-1
E1 (FD) : 3226.5732 cm^-1
Fundamental (FD) : 2144.8757 cm^-1
First 5 eigenvalues (cm^-1):
E0 = 1081.6974 cm^-1
E1 = 3226.5732 cm^-1
E2 = 5348.5220 cm^-1
E3 = 7449.8673 cm^-1
E4 = 9533.1324 cm^-1
======================================================================
COMPREHENSIVE COMPARISON OF ALL METHODS
======================================================================
Ground state energy (E0):
Harmonic : 1084.8131 cm^-1
PT1 : 1090.4999 cm^-1
PT2 : 1081.4940 cm^-1
Variational (N=20) : 1528.0451 cm^-1
Finite Difference : 1081.6974 cm^-1
Morse (reference) : 1081.5635 cm^-1
First excited state energy (E1):
Harmonic : 3254.4394 cm^-1
PT1 : 3282.8734 cm^-1
PT2 : 3224.6464 cm^-1
Variational (N=20) : 4547.9969 cm^-1
Finite Difference : 3226.5732 cm^-1
Morse (reference) : 3225.1929 cm^-1
Fundamental transition:
Harmonic : 2169.6262 cm^-1
PT1 : 2192.3734 cm^-1
PT2 : 2143.1524 cm^-1
Variational (N=20) : 3019.9518 cm^-1
Finite Difference : 2144.8757 cm^-1
Morse (reference) : 2143.6294 cm^-1
Error in fundamental vs Morse:
Harmonic : +25.9968 cm^-1
PT1 : +48.7440 cm^-1
PT2 : -0.4770 cm^-1
Variational (N=20) : +876.3224 cm^-1
Finite Difference : +1.2463 cm^-1
======================================================================
/var/folders/x2/ggtknn9s49s3dcq96ckcgfgh0000gn/T/ipykernel_43430/1192904635.py:107: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.
norm = np.sqrt(np.trapz(eigenvectors[:, i]**2, r))
/var/folders/x2/ggtknn9s49s3dcq96ckcgfgh0000gn/T/ipykernel_43430/1733439826.py:133: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.
V_nm = np.trapz(integrand, r)
Wavefunction comparison plots generated!