Tutorials for new students in computational chemistry
Tutorial 16. Interactive -- Morse Potential and diatomic spectrum database
Learning objectives
Learn basics about maniputlating datasets with the python package, Pandas
Learn to fit Morse potential based on experimental vibrational spectrum data
Try the interactive python code to visualize the Morse potential of low-lying electronic states of diatomic molecules
Background
In quantum mechanics class, we learn about the Morse potential, which is a widely used approximation to the potential of diatomic molecules.
\[V(r)=D_e[1-e^{-a(r-r_e)}]^2\]
How can we relate it to real vibrational spectrum data of diatomic molecules? You might have seen some example from the textbook: given the vibrational spectrum constants $\omega_e$, $\omega_ex_e$ of a diatomic molecule, we can obtain $D_e$ and $a$ in the Morse potential formula, and plot the Morse potential of that molecule.
Can we do better than just one textbook example?
Yes, we can visualize the Morse potential and the vibrational energy levels of hundreds of diatomic molecules in seconds! That’s simple with the following recipe:
An open-access diatomic spectrum dataset
Some basic knowledge about dataset manipulation
Some python code to draw Morse potential curves
These ingredients become available to us only in this era of big data.
Try and learn
Run the following code to visualize Morse potentials of hundreds of diatomic molecules based on spectrum data obtained from diatomic database.
Instructions
First, click “Activate” to activate the code block.
Once you see the buttons to “run” at the bottom left corner of the code block, click “run” to run the code.
Please be patient. Starting the kernel can be slow sometimes.
You will see a plot for the Morse potential of OH (two low lying electronic states).
Select other diatomic molecules from the dropdown menu to visualize their Morse potentials.
###########################################
# First load the dataset
# The dataset is downloaded from: https://rios.mp.fhi.mpg.de/
# Citation: Liu, X., Truppe, S., Meijer, G. et al. The diatomic molecular spectroscopy database. J Cheminform 12, 31 (2020). https://doi.org/10.1186/s13321-020-00433-8
############################################
!wget--trust-server-nameshttps://rios.mp.fhi.mpg.de/export_table.php-Omoleculedata.csv%matplotlibwidgetimportipywidgetsaswidgetsimportpandasaspdimportnumpyasnpimportmatplotlib.pyplotaspltimportmatplotlibasmplfromscipy.specialimportfactorialfromscipy.specialimportgenlaguerre,gammadf=pd.read_csv('moleculedata.csv')# Remove the records where the vibrational spectrum constants have NaN
df=df.dropna(subset=['omega_ex_e (cm^{-1})','omega_e (cm^{-1})'])# Create the dropdown list of diatomic molecules available
molecule_options=widgets.Dropdown(options=df['Molecule'].unique(),value='OH',description='Molecule:',disabled=False,)###########################################################################
# Now define the functions to plot Morse potential.
# The code is adapted from: https://scipython.com/blog/the-morse-oscillator/
# Major modifications to the code:
# 1. All member variables of the Morse class are stored in atomic units (energy: Hartree, mass: m_e, length: bohr, hbar=1) to simplify calculations and avoid unit conversions
# 2. Reduced mass as input instead of mA and mB
# 3. Plots have eV used for energy unit
# 4. Plot different electronic states of the same molecule on the same figure, label term symbol of the electronic states along the potential
##############################################################################
# Let's convert cm-1 to Eh (Hartree)
wavenumber2Hartree=4.55633e-6mp=1836# Porton mass is 1836 mass of electron
angstrom2Bohr=1.88973hbar=1hartree2eV=27.2114bohr2angstrom=0.529177classMorse:"""A class representing the Morse oscillator model of a diatomic."""def__init__(self,mu,we,wexe,re,Te):"""Initialize the Morse model for a diatomic molecule.
Input arguments:
mu input reduced mass (mass unit: proton mass).
we, wexe are the Morse parameters (cm-1).
re is the equilibrium bond length (Angstrom).
Te is the electronic energy (cm-1).
class member variables:(everything in atomic unit)
mu, in electron mass (standard in atomic unit)
we, wexe in Hartree
re in bohr
De in hartree
D0 in hatree
"""# First of all, convert everything to atomic unit
# where hbar = 1, me=1, Eh=1
self.mu=mu*mpself.we,self.wexe=we*wavenumber2Hartree,wexe*wavenumber2Hartreeself.re=re*angstrom2Bohrself.Te=Te*wavenumber2Hartreeself.De=self.we**2*hbar/4/self.wexe# Morse parameters, a and lambda.
self.a=self.calc_a()self.ke=2*self.De*self.a**2self.lam=np.sqrt(2*self.mu*self.De)/self.a/hbar# Maximum vibrational quantum number.
self.vmax=int(np.floor(self.lam-0.5))self.make_rgrid()self.V=self.Vmorse(self.r)defmake_rgrid(self,n=1000,rmin=None,rmax=None,retstep=False):"""Make a suitable grid of internuclear separations."""self.rmin,self.rmax=rmin,rmaxifrminisNone:# minimum r where V(r)=De on repulsive edge
self.rmin=self.re-np.log(2)/self.aifrmaxisNone:# maximum r where V(r)=f.De
f=0.999self.rmax=self.re-np.log(1-f)/self.aself.r,self.dr=np.linspace(self.rmin,self.rmax,n,retstep=True)ifretstep:returnself.r,self.drreturnself.rdefcalc_a(self):"""Calculate the Morse parameter, a.
Returns the Morse parameter, a, in unit bohr^-1
"""returnself.we*np.sqrt(self.mu/2.0/self.De)defVmorse(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).
"""returnself.De*(1-np.exp(-self.a*(r-self.re)))**2defEmorse(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.5return(self.we*vphalf-self.wexe*vphalf**2)defcalc_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 hatree). rm < rp.
"""b=np.sqrt(E/self.De)return(self.re-np.log(1+b)/self.a,self.re-np.log(1-b)/self.a)defcalc_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.
"""ifrisNone:r=self.rxi=2*self.lam*np.exp(-self.a*(r-self.re))alpha=2*(self.lam-v)-1psi=(xi**(self.lam-v-0.5)*np.exp(-xi/2)*genlaguerre(v,alpha)(xi))psi*=psi_max/np.max(psi)returnpsidefcalc_psi_Normalized(self,v,xi):alpha=2*(self.lam-v)-1psi=(xi**(self.lam-v-0.5)*np.exp(-xi/2)*genlaguerre(v,alpha)(xi))Nv=np.sqrt(factorial(v)*(2*self.lam-2*v-1)/gamma(2*self.lam-v))returnNv*psidefplot_V(self,ax,**kwargs):"""Plot the Morse potential on Axes ax."""ax.plot(self.r*bohr2angstrom,(self.V+self.Te)*hartree2eV,**kwargs)defget_vmax(self):"""Return the maximum vibrational quantum number."""returnint(self.we/2/self.wexe-0.5)defdraw_Elines(self,vlist,ax,**kwargs):"""Draw lines on Axes ax representing the energy level(s) in vlist."""ifisinstance(vlist,int):vlist=[vlist]forvinvlist:E=self.Emorse(v)rm,rp=self.calc_turning_pts(E)ax.hlines((E+self.Te)*hartree2eV,rm*bohr2angstrom,rp*bohr2angstrom,**kwargs)deflabel_levels(self,vlist,ax):ifisinstance(vlist,int):vlist=[vlist]forvinvlist:E=self.Emorse(v)rm,rp=self.calc_turning_pts(E)ax.text(s=r'$v={}$'.format(v),x=rp*bohr2angstrom+0.6,y=(E+self.Te)*hartree2eV,va='center')deflabel_electronstate(self,state_str,ax):E=self.Emorse(self.vmax)ax.text(s=state_str,x=np.max(self.r)*bohr2angstrom,y=(E+self.Te)*hartree2eV-0.2,va='center')defplot_psi(self,vlist,ax,r_plot=None,scaling=1,**kwargs):"""Plot the Morse wavefunction(s) in vlist on Axes ax."""ifisinstance(vlist,int):vlist=[vlist]forvinvlist:E=self.Emorse(v)*hartree2eVifr_plotisNone:rm,rp=self.calc_turning_pts(E)x=self.r[self.r<rp*1.2]else:x=r_plotpsi=self.calc_psi(v,r=x,psi_max=self.we/2)psi_plot=psi*scaling+(self.Emorse(v)+self.Te)*hartree2eVax.plot(x*bohr2angstrom,psi_plot,**kwargs)##########################################
# Below is the interactive plotting part #
##########################################
fig=plt.figure()@widgets.interact(molname=molecule_options)defupdate(molname):fig.clf()mol=df[df['Molecule']==molname]states=mol['Electronic state'].valuesax=plt.gca()# Set the default color cycle
colors=['blue','red','green','cyan','orange','gray','purple','brown','teal','coral']i_color=0forstateinstates:mycolor=colors[i_color]i_color+=1we=mol[mol['Electronic state']==state]['omega_e (cm^{-1})'].values[0]wexe=mol[mol['Electronic state']==state]['omega_ex_e (cm^{-1})'].values[0]re=mol[mol['Electronic state']==state]['Re (\AA)'].values[0]Te=mol[mol['Electronic state']==state]['Te (cm^{-1})'].values[0]mu=mol[mol['Electronic state']==state]['Reduced mass'].values[0]X=Morse(mu,we,wexe,re,Te)X.make_rgrid()X.V=X.Vmorse(X.r)X.plot_V(ax,color=mycolor)X.draw_Elines(range(X.vmax),ax,linewidth=0.5,color=mycolor)X.label_electronstate(state,ax)#X.plot_psi([2, 14], ax, scaling=2, color=COLOUR1)
#X.label_levels([2, 14], ax)
ax.set_xlabel(r'$r\;/\mathrm{\AA}$')ax.set_ylabel(r'Energy (eV)')ax.spines['top'].set_visible(False)ax.spines['right'].set_visible(False)
Questions
Read the code and answer the questions below
How do we obtain the binding energy $D_e$ from the vibrational spectrum constants $\omega_e$ and $\omega_e x_e$?
When the spectrum data for multiple electronic states are available, which variable in the data set determines the relative energy of different states?
What is the typical energy difference between electronic states? What is the typical energy difference between vibrational states?
Do you understand the meaning of the Term Symbols that label different electronic states?
I made some modifications to make it follow the conventions of writing quantum chemistry programs.
All member variables of the Morse class are stored in atomic units (energy: Hartree, mass: m_e, length: bohr, hbar=1) to simplify calculations and avoid unit conversions
Reduced mass is used as input instead of $m_A$ and $m_B$.
Plots use the unit, eV, for energy.
We plot different electronic states of the same molecule on the same figure, label term symbol of the electronic states along with the potential.