3
$\begingroup$

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:

Figure of harmonic and anharmonic PESs

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.

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ You can draw plots in LaTex. Check this link and this link to now further. $\endgroup$
    – Hemanth Haridas
    Jul 24 at 6:01
  • $\begingroup$ That doesn't really help me $\endgroup$
    – isolated matrix
    Jul 24 at 11:08
  • $\begingroup$ 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. $\endgroup$
    – Nike Dattani
    Jul 26 at 6:19
  • $\begingroup$ @NikeDattani I would like any answer that works $\endgroup$
    – isolated matrix
    Jul 27 at 4:25

1 Answer 1

Reset to default
4
$\begingroup$

A script that works for arbitrary potentials

My script has produced the following figures for this paper.

Simple single-well potential:

enter image description here

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

enter image description here

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.

Distinction among vibrational levels that have been spectroscopically observed versus ones that have not been observed:

enter image description here

Inset showing vibrational levels:

enter image description here

Vibrational levels labeled:

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

enter image description here

It was also used for the figure in this paper.

Selected vibrational levels are shown or not shown:

enter image description here

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);
Improve this answer
$\endgroup$
5
  • 1
    $\begingroup$ 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 $\endgroup$
    – isolated matrix
    Jul 25 at 4:50
  • 2
    $\begingroup$ @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)$. $\endgroup$
    – Nike Dattani
    Jul 25 at 13:19
  • $\begingroup$ 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... $\endgroup$
    – isolated matrix
    Jul 26 at 5:15
  • $\begingroup$ 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. $\endgroup$
    – Nike Dattani
    Jul 26 at 6:24
  • $\begingroup$ 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. $\endgroup$
    – isolated matrix
    Jul 26 at 13:31

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .