Script to draw one-dimensional PES comparing harmonic and anharmonic vibrational modes

This is an extension of a question I asked a little while ago: Software to draw one-dimensional PES including vibrational energy levels. I am asking the followup here based on advice received after asking the other question.

Again, my aim is to compare the difference between a harmonic and anharmonic oscillator. While I am definitely able to produce a Morse potential using spectroscopic constants (two very good Python scripts can be found in the answers to the previous question), I'm not sure how I can generate an equivalent diagram for the harmonic PES.

I believed that it would actually be more simple than the Morse function (after all, the vibrational levels are equally spaced and the function is, essentially, a parabola). However, the only method of generating such a diagram is another Python script (The Harmonic Oscillator Wave Functions) which is, to my eyes, far more complicated than the script the same person provided for Morse PESs, to the extent that I cannot figure out where to insert the experimental spectroscopic parameters...

While looking for an answer, though, I came across this figure from LibreTexts:

This is exactly the type of figure I would like to generate, as it clearly shows the difference between the harmonic and anharmonic vibrational energy levels. It should be possible to plot both of the PESs given the same spectroscopic constants for the same molecule, but I can't figure out how to manipulate the scripts to do so.

So, does anyone know of a way of generating a harmonic PES with vibrational energy levels marked, similar to the answers given for the Morse potential (Software to draw one-dimensional PES including vibrational energy levels), or, even better, does anyone know how I could draw a diagram like the one in the figure above? I have the anharmonicity and force constants ($$\omega$$, $$\omega x$$, etc.), but don't have a list of frequencies.

• You can draw plots in LaTex. Check this link and this link to now further. Commented Jul 24, 2023 at 6:01
• That doesn't really help me Commented Jul 24, 2023 at 11:08
• If you want to use/modify that specific script that you've mentioned in your edit well after my answer, I recommend to ask a new question. Commented Jul 26, 2023 at 6:19
• @NikeDattani I would like any answer that works Commented Jul 27, 2023 at 4:25

A script that works for arbitrary potentials

My script has produced the following figures for this paper.

Double-minimum potential with a small barrier and a "shelf":

Notice that when the potential gets wider, the vibrational levels get closer together, as we expect based on what happens to a "particle in a box" when the box gets wider.

The same basic script was used for this paper.

Vibrational levels labeled:

The same basic script was also used for the diagram in my MMSE answer here.

It was also used for the figure in this paper.

Here is the script that was used for the first paper listed above

close('all')
figure1=figure(1);

set(gca,'Position',[0.101190476190476 0.105359342915811 0.885 0.888357289527721])
set(gcf,'Color','w')

minX=1;maxX=20;
splineMeshSize_for_r=0.00001; % I just chose it to be of the same precision as r_e defined in cell 1. Splines can only be necessary for short-range, since long-range seems to have points very close together.
splineMesh_for_r=minX:splineMeshSize_for_r:maxX;

V_ab_initio_spline_6X=spline(r_abInitio_6X,V_abInitio_6X,splineMesh_for_r);
V_ab_initio_spline_6a=spline(r_abInitio_6a,V_abInitio_6a,splineMesh_for_r);
V_ab_initio_spline_3A=spline(r_abInitio_3A,V_abInitio_3A,splineMesh_for_r);
V_ab_initio_spline_3B=spline(r_abInitio_3B,V_abInitio_3B,splineMesh_for_r);
V_ab_initio_spline_3C=spline(r_abInitio_3C,V_abInitio_3C,splineMesh_for_r);
V_ab_initio_spline_3b=spline(r_abInitio_3b,V_abInitio_3b,splineMesh_for_r);
V_ab_initio_spline_3c=spline(r_abInitio_3c,V_abInitio_3c,splineMesh_for_r);
V_ab_initio_spline_3d=spline(r_abInitio_3d,V_abInitio_3d,splineMesh_for_r);

minY=-6000;maxY=300;
axis([minX,maxX,minY,maxY])
line([minX maxX], [0 0],'Color','k','LineWidth',3)

hold('on')

% [Vmin_ab_initio, indexOfVmin_ab_initio]=min(V_ab_initio_spline);
% plotPointsSymmetricallyOnPotential(vibrationalEnergies(1:2:end),V_ab_initio_spline,splineMesh_for_r,indexOfVmin_ab_initio);

axis([minX,maxX,minY,maxY])
plot(r_abInitio_6X,V_abInitio_6X,'Color',[0,255,100]./255,'Linewidth',3);
plot(r_abInitio_3A,V_abInitio_3A,'Color','r','Linewidth',3);
plot(r_abInitio_3B,V_abInitio_3B,'Color','k','Linewidth',3);
plot(r_abInitio_3C,V_abInitio_3C,'Color','c','Linewidth',3);

figure2=figure(2);

set(gca,'Position',[0.101190476190476 0.105359342915811 0.885 0.888357289527721])
set(gcf,'Color','w')

minX=1;maxX=20;
minY=-10000;maxY=300;
axis([minX,maxX,minY,maxY])
line([minX maxX], [0 0],'Color','k','LineWidth',3)

splineMeshSize_for_r=0.00001; % I just chose it to be of the same precision as r_e defined in cell 1. Splines can only be necessary for short-range, since long-range seems to have points very close together.
splineMesh_for_r=minX:splineMeshSize_for_r:maxX;

hold('on')

plot(r_abInitio_6a,V_abInitio_6a,'Color','b','Linewidth',3);
plot(r_abInitio_3b,V_abInitio_3b,'Color','m','Linewidth',3);
plot(r_abInitio_3c,V_abInitio_3c,'Color','g','Linewidth',3);
plot(r_abInitio_3d,V_abInitio_3d,'Color',[100,100,100]./255,'Linewidth',3);

% [Vmin_ab_initio, indexOfVmin_ab_initio]=min(V_ab_initio_spline);
% plotPointsSymmetricallyOnPotential(vibrationalEnergies(1:2:end),V_ab_initio_spline,splineMesh_for_r,indexOfVmin_ab_initio);

axis([minX,maxX,minY,maxY])
plot(r_abInitio_6a,V_abInitio_6a,'Color','b','Linewidth',3);
plot(r_abInitio_3b,V_abInitio_3b,'Color','m','Linewidth',3);
plot(r_abInitio_3c,V_abInitio_3c,'Color','g','Linewidth',3);
plot(r_abInitio_3d,V_abInitio_3d,'Color',[100,100,100]./255,'Linewidth',3);

%%
close('all');figure3=figure(3);hold('on')
set(gca,'Position',[0.101190476190476 0.114161849710983 0.892953008436171 0.879554782732549])
set(gcf,'Color','w')

minX=2;maxX=15;
minY=-6000;maxY=300;
axis([minX,maxX,minY,maxY])
line([minX maxX], [0 0],'Color','k','LineWidth',3)
re=3.98;

plot(r_abInitio_3b,V_abInitio_3b,'Color',[0,255,100]./255,'Linewidth',10);
plot(r_3b_p6q7r05_90,V_3b_p6q7r05_90,'k','LineWidth',3)

v0=-5654.932449854602055;
v( 1)=-5493.765135128664951;
v( 2)=-5385.710984514686970;
v( 3)=-5241.269953247034209;
v( 4)=-5083.976206177168933;
v( 5)=-4914.85;
v( 6)=-4737.691247106215087;
v( 7)=-4554.829936403927604;
v( 8)=-4368.240502184650722;
v( 9)=-4179.292096277415112;
v(10)=-3988.756231213321371;
v(11)=-3796.972235927646125;
v(12)=-3604.218750080568952;
v(13)=-3411.033149015913750;
v(14)=-3218.297919388399805;
v(15)=-3027.108781656958399;
v(16)=-2838.477232599816034;
v(17)=-2652.930601841349471;
v(18)=-2470.302596380747673;
v(19)=-2290.105676375571875;
v(20)=-2112.256988503596288;
v(21)=-1937.405398146945536;
v(22)=-1766.631273238948864;
v(23)=-1600.818435795562496;
v(24)=-1440.124051819757312;
v(25)=-1284.340696145073408;
v(26)=-1134.583144744511360;
v(27)= -996.778888028432000;
v(28)= -914.777113202062848;
v(29)= -882.576799450836608;
v(30)= -845.026425674076544;
v(31)= -800.327833982061696;
v(32)= -753.924510206392448;
v(33)= -705.900692148498688;
v(34)= -654.935744351130112;
v(35)= -602.077023609842304;
v(36)= -548.234951970257472;
v(37)= -493.706925401166272;
v(38)= -438.926074728071424;
v(39)= -384.474432243322304;
v(40)= -330.903787772794240;
v(41)= -278.814572462308608;
v(42)= -228.913912424278848;
v(43)= -182.018345069401152;
v(44)= -139.060176467093648;
v(45)= -101.066168262798576;
v(46)=  -69.047731355521066;
v(47)=  -43.749645705548336;
v(48)=  -25.299916111188852;
v(49)=  -13.022995406988458;
v(50)=   -5.669747141689355;
v(51)=   -1.855271158639904;
v(52)=   -0.324750667836605;
v(53)=   -0.002393731491646;
vibrationalEnergies=[v0 v];

[Vmin_ab_initio_3b, indexOfVmin_ab_initio_3b]=min(V_ab_initio_spline_3b);
%plotPointsSymmetricallyOnPotential(vibrationalEnergies(1:2:end),V_ab_initio_spline_3b,splineMesh_for_r,indexOfVmin_ab_initio_3b,200);

vibrationalLevelsMeasured=[0:53];
for vibrationalLevel=vibrationalLevelsMeasured+1
line([interp1(V_3b_p6q7r05_90(r_3b_p6q7r05_90<re),r_3b_p6q7r05_90(r_3b_p6q7r05_90<re),vibrationalEnergies(vibrationalLevel)) interp1(V_3b_p6q7r05_90(r_3b_p6q7r05_90>re),r_3b_p6q7r05_90(r_3b_p6q7r05_90>re),vibrationalEnergies(vibrationalLevel))], [vibrationalEnergies(vibrationalLevel) vibrationalEnergies(vibrationalLevel)],'Color','b','LineWidth',3)
r_innerTurningPoint(vibrationalLevel)=interp1(V_3b_p6q7r05_90(r_3b_p6q7r05_90<re),r_3b_p6q7r05_90(r_3b_p6q7r05_90<re),vibrationalEnergies(vibrationalLevel));
scatter(r_innerTurningPoint(vibrationalLevel),vibrationalEnergies(vibrationalLevel),'MarkerFaceColor','b','MarkerEdgeColor','b');
r_outerTurningPoint(vibrationalLevel)=interp1(V_3b_p6q7r05_90(r_3b_p6q7r05_90>re),r_3b_p6q7r05_90(r_3b_p6q7r05_90>re),vibrationalEnergies(vibrationalLevel));
scatter(r_outerTurningPoint(vibrationalLevel),vibrationalEnergies(vibrationalLevel),'MarkerFaceColor','b','MarkerEdgeColor','b');
end

plot(r(700:end),V(end,700:end),'r','LineWidth',6)

annotation(figure3,'textbox',[0.620 0.760 0.327 0.087],'Color','r','String','$$V(r) = - \frac{D_3(r)C_3}{r^3} - \frac{D_6(r)C_6}{r^6} - \frac{D_8C_8(r)}{r^8}$$','LineStyle','none','Interpreter','latex','FontSize',24);
annotation(figure3,'textbox',[0.421 0.609 0.107 0.0636],'String','MLR$$_{5.9}^{6,7}(17)$$','LineStyle','none','Interpreter','latex','FontSize',24);
annotation(figure3,'textbox',[0.293 0.223 0.107 0.0636],'String','Li$$_2\left(3b,3^3\Pi_u\right)$$','LineStyle','none','Interpreter','latex','FontSize',36,'FontName','Helvetica','FitBoxToText','off');
annotation(figure3,'textbox',[0.438 0.532 0.107 0.0636],'Color',[0,255,100]./255,'String','Original','LineStyle','none','Interpreter','latex','FontSize',24);
annotation(figure3,'line',[0.366764275256223 0.412884333821376],[0.63775968992248 0.637209302325581],'LineWidth',3);
scatter(5.89,-2750,200,'MarkerFaceColor',[0,255,100]./255,'MarkerEdgeColor','k')
scatter(6.5,-2750,200,'MarkerFaceColor',[0,255,100]./255,'MarkerEdgeColor','k')
annotation(figure3,'arrow',[0.620790629575403 0.588579795021962],[0.850162790697674 0.924031007751938],'LineWidth',3);
annotation(figure3,'arrow',[0.922401171303075 0.953147877013177],[0.768992248062015 0.689922480620155],'LineWidth',3);

xlabel('Internuclear distance \AA','Interpreter','Latex','FontSize',36)                                                                                    % (r) taken out as per request by Kirk Madison in email on 25/8/2012
ylabel('$$V(r)$$ cm$$^{-1}$$','Interpreter','Latex','FontSize',36)
box('on');
set(gca,'XMinorTick','on','YMinorTick','on','LineWidth',2,'FontSize',20);

• This looks amazing! However... do you have to have a list of the frequencies to use this script? I don't have that – I only have the anharmonicity constants (I forgot to add that detail to the question... I'll do so now). Also, my question was more about generating a harmonic PES, not an anharmonic PES Commented Jul 25, 2023 at 4:50
• @isolatedmatrix the script doesn't solve the Schrödinger equation, so the energy levels need to be provided from somewhere, however for a harmonic potential the energy levels are just $\hbar\omega(n+1/2)$. Commented Jul 25, 2023 at 13:19
• I know that, but I don't know how to generate the figure. I can absolutely calculate the vibrational levels if I need to, but what I can't do is draw the parabola with the vibrational levels marked, like in your answer for the anharmonic potential surface. The harmonic oscillator script I linked above can do it, apparently, but I have no idea how to modify it to a) centre the parabola on the equilibrium bond length; b) remove the wave functions; c) customise it to apply to the O$_2$ molecule's vibrational levels. I might add that to the question, actually... Commented Jul 26, 2023 at 5:15
• In my script, you can make V(r) any parabola you want. However if you want to work with a specific pre-existing script, please post a new question, this time specifically aimed at modifying that script. Commented Jul 26, 2023 at 6:24
• If I can make V(r) a parabola, that would be good, but I don't know how to modify your script to do that (I don't understand Python). I honestly don't even know where the command to plot anything even is. Commented Jul 26, 2023 at 13:31