I'm sorry if this is a really super basic question. I am looking at the potential curves of a diatomic molecule in a book , and to plot them. I google the molecule and got to the NIST website, there I get diatomic mass and constants for each state, such as:

  • Te - minimum electronic energy (cm-1)
  • ωe - vibrational constant – first term (cm-1)
  • ωexe - vibrational constant – second term (cm-1)
  • ωeye - vibrational constant – third term (cm-1)
  • Be - rotational constant in equilibrium position (cm-1)
  • αe - rotational constant – first term (cm-1)
  • γe - rotation-vibration interaction constant (cm-1)
  • De - centrifugal distortion constant (cm-1)
  • βe - rotational constant – first term, centrifugal force (cm-1)
  • re - internuclear distance (Å)
  • Trans. - observed transition(s) corresponding to electronic state

Great, I should be able to use these to plot the potential energy curves for each electronic state somehow right? there is no formula i can see that here, is it just morse potential? what happens if the state is unbound? Please help me understand how to plot the potential curves. Thanks!

  • 1
    $\begingroup$ You could try the method in this paper. sciencedirect.com/science/article/abs/pii/009784859500016L $\endgroup$
    – porphyrin
    Commented Oct 3, 2020 at 7:45
  • $\begingroup$ Im sorry, that doesn't answer the question. and that paper is from 1995. mathematica (which I dont have) is probably very different today. the paper related to RKR potentials, why is that appropriate for for this question? $\endgroup$
    – dpdp
    Commented Oct 4, 2020 at 7:11
  • 1
    $\begingroup$ The fact that it is Mathematica is unimportant. The equations for the classical turning points of the potential for the potential are given in terms of the parameters you quote. I coded this up in Python (free) quite easily. If you cannot see equations due to pay wall I will put them in an answer $\endgroup$
    – porphyrin
    Commented Oct 4, 2020 at 7:21
  • $\begingroup$ yes, I can't access that paper. Thanks for the clarification. I'd very appreciate if you can help with posting the equation in your answer. $\endgroup$
    – dpdp
    Commented Oct 5, 2020 at 17:19

1 Answer 1


The RKR method determines the classical turning points $r_{\pm}$ of the potential energy based on knowing experimentally determined spectroscopic constants. ( Rydberg, 1931, Z. Physik. v73,376, Klein, 1932 Z. Physik. v76, 226, Rees Proc. Phys. Soc.(Lond.) 1947, v59, 998.)

The energy equation takes the usual form. Parameters all in cm$^{-1}$, distances in m, reduced mass $\mu$ in u, $h$ in cm$^{-1}$s.

$$E_v=\omega_e(v+1/2) - x_e\omega_e(v+1/2)^2 + y_e\omega_e(v+1/2)^3 + z_e\omega_e(v+1/2)^4 +\cdots$$

The vibrational quantum number $v$ should be considered as a variable that has values from $-1/2$ to whatever positive value is required.

The derivative of the rotational energy wrt. $J$ at $J=0$ is also needed,

$$B'(v)= B_e-\alpha_e(v+1/2)+\gamma_e(v+1/2)^2+\cdots$$

The function representing the turning points is

$$r_\pm= \frac{f(v)}{2}\left (\sqrt{1+\frac{4}{f(v)g(v)}}\pm 1 \right)$$

where functions $f, \, g$ are integrals. These functions 'blow up' at the upper integration limit $x \to v$ so a correction is added to prevent this. The basic equation is

$$f(v)=\frac{h}{\pi\sqrt{2\mu}}\int_{-1/2}^v\frac{dx}{ \sqrt{E(v)-E(x)}}$$

and the modified one

$$f(v)=\frac{h}{\pi\sqrt{2\mu}}\left(\int_{-1/2}^{v-\delta}\frac{dx}{ \sqrt{E(v)-E(x)} }+ 2\sqrt{\frac{\delta}{Q_v}}\right)$$

where $\delta$ is small $\approx 10^{-5}$ and

$$Q_v=\omega_e - 2x_e\omega_e(v+1/2) + 3y_e\omega_e(v+1/2)^2 + 4z_e\omega_e(v+1/2)^3 +\cdots$$

The function $g$ with correction is

$$g(v)=\frac{4\pi\sqrt{2\mu}}{h} \left( \int_{-1/2}^{v-\delta} \frac{B'(x)}{\sqrt{E(v)-E(x)} }dx+2B'(v)\sqrt{\frac{\delta}{Q_v} } \right) $$

A set of values for CO is
$\mu = 6.85841,\; \omega_e= 2169.8135,\; x_e\omega_e=13.2883,\; y_e\omega_e=0.010511,\; z_e\omega_e= 0.000057$ $ B_e= 1.9312808,\; \alpha_e= 0.0175044,\; \gamma_e= 0.000000548$.

Plot energy $E(x)$(cm$^{-1}) $ vs $r_\pm(x)$ (pm) is shown. The value at $x=-1/2$ should be $r_e=103.14$ pm the equilibrium bond length with energy of zero. As values are based on cm$^{-1}$ a scaling of $10^{14}$ is needed to go to pm.


I used Jupyter notebooks (via Anaconda) and python 3 to do the calculation using built in integrator from numpy/scipy and matplotlib to draw the graph. All are free and very easy to use.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.