Implementation

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import assoc_laguerre, gamma

class MorseOscillator:

    def __init__(self, m, a, De, re):
        self.mass = m
        self.a = a
        self.De = De
        self.re = re
        
        self.lam = np.sqrt(2.0*self.mass*self.De)/self.a
        self.omega0 = self.a * np.sqrt(2.0*self.De/self.mass)
    
    def get_potential(self, r):
        return self.De*(1.0 - np.exp(-self.a*(r - self.re)))**2

    def get_energy(self, n):
        return self.omega0*(n + 0.5) - (self.omega0*(n + 0.5))**2/(4.0*self.De)

    def get_eigenfunction(self, r, n):
        z = 2.0*self.lam*np.exp(-(r - self.re))
        norm = np.sqrt((np.math.factorial(n)*(2.0*self.lam - 2.0*n - 1)) \
            / (gamma(2.0*self.lam - n)))
        p = assoc_laguerre(z, n, k=(2.0*self.lam - 2.0*n - 1))
        return norm * z**(self.lam - n - 0.5) * np.exp(-0.5*z) * p
    
    def plot(self, xvalues=np.linspace(0.5, 8, 1000), nstate=17, scale=.01):
        fig, ax = plt.subplots(figsize=(10, 6))
        ax.plot(xvalues, self.get_potential(xvalues), lw=5.0, color="blue")
        plt.xlabel(r"$x$")
        plt.ylabel(r"$V(x)$")
        for state in range(nstate):
            en = self.get_energy(state)
            ax.plot([xvalues[0], xvalues[-1]], [en,en], lw=2.0, color="red")
            ax.plot(xvalues, en + scale * self.get_eigenfunction(xvalues, state),
                    lw=2.0, color="black", alpha=.5)
        return fig, ax

As a test case, we consider Morse potential for the \(\mathrm{H_{2}}\) molecule that has the following parameters (in atomic units): $$ \begin{align} D_{e} &= 0.1744 \\ a &= 1.02764 \\ x_{e} &= 1.40201 \end{align} $$

The mass is defined as the reduced mass of a diatomic molecule, which, in the case of \(\mathrm{H_{2}}\), corresponds to the half of the atomic mass of hydrogen,

$$ m = \frac{m_{H}}{2}.$$

me = 9.1093837015e-31
u = 1.66053906660e-27
auofmass = u / me

De = 0.1744
a = 1.02764
xe = 1.40201
m = (1.0078250322 / 2.) * auofmass

h2morse = MorseOscillator(m, a, De, xe)
fig, ax = h2morse.plot()
ax.set_ylim(-.01, 0.19)
ax.grid()

plt.show()