The Morse oscillator


The Morse oscillator is a model for a vibrating diatomic molecule that improves on the simple harmonic oscillator model in that the vibrational levels converge with increasing energy and that at some finite energy the molecule dissociates. The potential energy varies with displacement of the internuclear separation from equilibrium, $x = r - r_\mathrm{e}$ as: $$ V(x) = D_\mathrm{e}\left[ 1-e^{-ax} \right]^2, $$ where $D_\mathrm{e}$ is the dissociation energy, $a = \sqrt{k_\mathrm{e}/2D_\mathrm{e}}$, and $k_\mathrm{e} = (\mathrm{d}^2V/\mathrm{d}x^2)_\mathrm{e}$ is the bond force constant at the bottom of the potential well.

enter image description here

The Morse oscillator Schrödinger equation, $$ -\frac{\hbar^2}{2m}\frac{\mathrm{d}^2\psi}{\mathrm{d}x^2} + V(x)\psi = E\psi $$ can be solved exactly.

It is helpful to define the new parameters, $$ \lambda = \frac{\sqrt{2mD_\mathrm{e}}}{a\hbar}\quad \mathrm{and} \quad z = 2\lambda e^{-x}, $$ in terms of which the eigenfunctions are $$ \psi_v(z) = N_v z^{\lambda - v - \frac{1}{2}}e^{-z/2}L_v^{(2\lambda - 2v - 1)}(z). $$ Here, $$ L_v^{(\alpha)} = \frac{z^{-\alpha}e^z}{v!}\frac{\mathrm{d}^v}{\mathrm{d}z^v}(z^{v+\alpha}e^{-z}) $$ is a generalized Laguerre polynomial and $$ N_v = \sqrt{\frac{v!(2\lambda - 2v - 1)}{\Gamma(2\lambda-v)}} $$ is a normalization constant. The quantum number, $v$ takes values $$ v = 0, 1, 2, \cdots, \lfloor \lambda - \frac{1}{2} \rfloor. $$ The corresponding eigenvalues are usually written $$ E_v = \omega_\mathrm{e}(v+\frac{1}{2}) - \omega_\mathrm{e}x_\mathrm{e}(v+\frac{1}{2})^2, $$ where $$ \omega_\mathrm{e} = \frac{a}{2\pi c}\sqrt{\frac{2D_\mathrm{e}}{m}} \quad \mathrm{and} \quad \omega_\mathrm{e} x_\mathrm{e} = \frac{\omega_\mathrm{e} ^2}{4D_\mathrm{e} }. $$

The code below defines the class Morse representing a Morse oscillator, which can be used to generate depictions of the wavefunctions and energy levels as given in the example for $\mathrm{^1H^{35}Cl}$ above.

import numpy as np
from scipy.constants import h, hbar, c, u
from scipy.misc import factorial
from scipy.special import genlaguerre, gamma

# Factor for conversion from cm-1 to J
FAC = 100 * h * c

class Morse:
    """A class representing the Morse oscillator model of a diatomic."""

    def __init__(self, mA, mB, we, wexe, re, Te):
        """Initialize the Morse model for a diatomic molecule.

        mA, mB are the atom masses (atomic mass units).
        we, wexe are the Morse parameters (cm-1).
        re is the equilibrium bond length (m).
        Te is the electronic energy (minimum of the potential well; origin
            of the vibrational state energies).


        self.mA, self.mB = mA, mB = mA*mB/(mA+mB) * u
        self.we, self.wexe = we, wexe = re
        self.Te = Te

        self.De = we**2 / 4 / wexe * FAC = (2 * np.pi * c * 100 * we)**2 *
        #  Morse parameters, a and lambda.
        self.a = self.calc_a()
        self.lam = np.sqrt(2 * * self.De) / self.a / hbar
        # Maximum vibrational quantum number.
        self.vmax = int(np.floor(self.lam - 0.5))

        self.V = self.Vmorse(self.r)

    def make_rgrid(self, n=1000, rmin=None, rmax=None, retstep=False):
        """Make a suitable grid of internuclear separations."""

        self.rmin, self.rmax = rmin, rmax
        if rmin is None:
            # minimum r where V(r)=De on repulsive edge
            self.rmin = - np.log(2) / self.a
        if rmax is None:
            # maximum r where V(r)=f.De
            f = 0.999
            self.rmax = - np.log(1-f)/self.a
        self.r, self.dr = np.linspace(self.rmin, self.rmax, n,
        if retstep:
            return self.r, self.dr
        return self.r

    def calc_a(self):
        """Calculate the Morse parameter, a.

        Returns the Morse parameter, a, from the equilibrium
        vibrational wavenumber, we in cm-1, and the dissociation
        energy, De in J.


        return (self.we * np.sqrt(2 * * np.pi *
                c * 100)

    def Vmorse(self, r):
        """Calculate the Morse potential, V(r).

        Returns the Morse potential at r (in m) for parameters De
        (in J), a (in m-1) and re (in m).


        return self.De * (1 - np.exp(-self.a*(r -**2

    def Emorse(self, v):
        """Calculate the energy of a Morse oscillator in state v.

        Returns the energy of a Morse oscillator parameterized by
        equilibrium vibrational frequency we and anharmonicity
        constant, wexe (both in cm-1).

        vphalf = v + 0.5
        return (self.we * vphalf - self.wexe * vphalf**2) * FAC

    def calc_turning_pts(self, E):
        """Calculate the classical turning points at energy E.

        Returns rm and rp, the classical turning points of the Morse
        oscillator at energy E (provided in J). rm < rp.


        b = np.sqrt(E / self.De)
        return ( - np.log(1+b) / self.a,
       - np.log(1-b) / self.a)

    def calc_psi(self, v, r=None, normed=True, psi_max=1):
        """Calculates the Morse oscillator wavefunction, psi_v.

        Returns the Morse oscillator wavefunction at vibrational
        quantum number v. The returned function is "normalized" to
        give peak value psi_max.


        if r is None:
            r = self.r
        z = 2 * self.lam * np.exp(-self.a*(r -
        alpha = 2*(self.lam - v) - 1
        psi = (z**(self.lam-v-0.5) * np.exp(-z/2) *
               genlaguerre(v, alpha)(z))
        psi *= psi_max / np.max(psi)
        return psi

    def calc_psi_z(self, v, z):
        alpha = 2*(self.lam - v) - 1
        psi = (z**(self.lam-v-0.5) * np.exp(-z/2) *
               genlaguerre(v, alpha)(z))
        Nv = np.sqrt(factorial(v) * (2*self.lam - 2*v - 1) /
                     gamma(2*self.lam - v))
        return Nv * psi

    def plot_V(self, ax, **kwargs):
        """Plot the Morse potential on Axes ax."""

        ax.plot(self.r*1.e10, self.V / FAC + self.Te, **kwargs)

    def get_vmax(self):
        """Return the maximum vibrational quantum number."""

        return int(self.we / 2 / self.wexe - 0.5)

    def draw_Elines(self, vlist, ax, **kwargs):
        """Draw lines on Axes ax representing the energy level(s) in vlist."""

        if isinstance(vlist, int):
            vlist = [vlist]
        for v in vlist:
            E = self.Emorse(v)
            rm, rp = self.calc_turning_pts(E)
            ax.hlines(E / FAC + self.Te, rm*1.e10, rp*1e10, **kwargs)

    def label_levels(self, vlist, ax):
        if isinstance(vlist, int):
            vlist = [vlist]

        for v in vlist:
            E = self.Emorse(v)
            rm, rp = self.calc_turning_pts(E)
            ax.text(s=r'$v={}$'.format(v), x=rp*1e10 + 0.6,
                    y=E / FAC + self.Te, va='center')

    def plot_psi(self, vlist, ax, r_plot=None, scaling=1, **kwargs):
        """Plot the Morse wavefunction(s) in vlist on Axes ax."""
        if isinstance(vlist, int):
            vlist = [vlist]
        for v in vlist:
            E = self.Emorse(v)
            if r_plot is None:
                rm, rp = self.calc_turning_pts(E)
                x = self.r[self.r<rp*1.2]
                x = r_plot
            psi = self.calc_psi(v, r=x, psi_max=self.we/2)
            psi_plot = psi*scaling + self.Emorse(v)/FAC + self.Te
            ax.plot(x*1.e10, psi_plot, **kwargs)

The following code generates the plot above for $\mathrm{^1H^{35}Cl}$:

import numpy as np
from matplotlib import rc
import matplotlib.pyplot as plt
from scipy.constants import h, c
from morse import Morse, FAC
rc('font', **{'family': 'serif', 'serif': ['Computer Modern'], 'size': 14})
rc('text', usetex=True)

COLOUR1 = (0.6196, 0.0039, 0.2588, 1.0)

# Atom masses and equilibrium bond length for (1H)(35Cl).
mA, mB = 1., 35.
X_re = 1.27455e-10
X_Te = 0
X_we, X_wexe = 2990.945, 52.818595

X = Morse(mA, mB, X_we, X_wexe, X_re, X_Te)
X.V = X.Vmorse(X.r)

fig, ax = plt.subplots()
X.plot_V(ax, color='k')

X.draw_Elines(range(X.vmax), ax)
X.draw_Elines(X.get_vmax(), ax, linestyles='--', linewidths=1)
X.plot_psi([2, 14], ax, scaling=2, color=COLOUR1)
X.label_levels([2, 14], ax)


Current rating: 4.2


Comments are pre-moderated. Please be patient and your comment will appear soon.

Hugh Jazz 3 years, 9 months ago

That's brilliant ! Do you know how to manage to plot two potentials of different states into one diagram ? I need to visualize vibronic transitions between two states of molecular nitrogen.

Link | Reply
Current rating: 5

christian 3 years, 9 months ago

I'm glad you like it!

You should be able to create two Morse objects (with different parameters, and Te being their electronic energy separation) and then plot them to the same Axes object, ax, by passing ax to plot_V, plot_psi, etc...

I might try to work up an example if you like.

Link | Reply
Currently unrated

Hugh Jazz 3 years, 9 months ago

Thanks for your reply ! It would brighten up my day if you could give me an example !

Link | Reply
Currently unrated

christian 3 years, 9 months ago

Check out

Link | Reply
Currently unrated

Hugh Jazz 3 years, 9 months ago

This is great ! I highly appreciate your work and I'm definitely going to purchase your book !

Link | Reply
Current rating: 5

christian 3 years, 9 months ago

A great choice! :)

Link | Reply
Current rating: 5

New Comment


required (not published)