Quantum Theoretical Chemistry, NWI-MOL112 (3EC)
Computer assignment 3: variational calculation atom-diatom system
Read the entire assignment before you start coding (29-Mar-2020).
Compute the bending modes of the Ar-CO complex with a variational
calculation using a coupled angular momentum basis \( |(jl)JM> \).
We will do \(J=0, M=0\).
Use the following data:
- Fix the CO distance at 1.128*Ang
Ang = 1.8897261 % bohr
- Fix the distance between Ar and the center-of-mass of CO at 3.56*Ang
- The bending potential (in hartree), expanded in Legendre polynomials,
is given by
V_bend(z) = sum_{L=0}^6 c_L P_L(z)
z = cos(theta)
cL=[ -6.77755480e-06 % c_0
1.46780154e-03 % c_1
1.69566457e-03 % c_2
9.44441298e-04 % c_3
6.24345128e-04 % c_4
2.32520826e-04 % c_5
1.05993802e-04] % c_6
The mases in amu are:
amu = 1822.8885 % m_e
mC = 12 * amu
mO = 15.994915 * amu
mAr = 39.962383 * amu
- A Matlab function to compute Clebsch-Gordan coefficients
is available here:
clebsch_gordan.m (bug fixed on
March 29, 2020)
Hints and intermediate results to help debug your code
Plot potential and find minimum
The bending potential is a linear combination of Legendre polynomials
\[
V_\mathrm{\,bend}(z) = \sum_{L=0}^6 c_L P_L(z).
\]
In Matlab, the Legendre polynomials can be computed with "legendreP(Ls, z)",
where Ls can be an array of L's. You can also use the recursion relations
from Chapter 8.1 of the lecture notes. These are some points to test
your potential in cm\(^{-1}\)
z V_bend(z)
_____________________
-1.00000 -49.49981
-0.50000 -155.06913
0.00000 -143.44938
0.50000 -5.12484
1.00000 1111.41725
Check that the minimum of the potential is \(V_\mathrm{min} = -159.34\) cm\(-1\).
Matrix elements in the uncoupled basis
Before you compute the matrix elements in the coupled basis, I
recommend writing a function that can compute matrix elements
of Racah spherical harmonics in a basis of Spherical harmonics, i.e.,
implement Eq. (8.7) from the lecture notes:
\[
\langle j_1 m_1 | C_{LM_L} | j_2 m_2 \rangle =
\sqrt{\frac{2j_2 + 1}{2j_1+1}}\;\langle L M_L j_2 m_2| j_1 m_1\rangle
\langle L, 0, j_2, 0| j_1, 0\rangle.
\]
Here are some results to test your function:
\[
\langle 2, 0|C_{3,0}|1 0\rangle = 0.3319700011
\]
\[
\langle 2, 2|C_{3,1}|1 1\rangle = -0.1106566670
\]
Potential energy matrix elements in the coupled basis
In order to compute the matrix elements in the coupled basis,
the potential must be written in the ``space-fixed'' coordinate
system using the Spherical Harmonics Addition Theorem,
see question 1c of week 5. I recommend first testing the
function to compute Clebsch-Gordan coefficients, e.g.,
by comparing to the
table on Wikipedia.
For total angular momentum \(J=0\), we must have \(j=l\) in the basis. In
a basis with \(j=0,1,2,3\) (and \(J=0\)), the potential energy matrix in
cm\(^-1\) is
V =
-1.4875 -185.9906 166.4329 -78.3448
-185.9906 147.3746 -235.1661 186.0490
166.4329 -235.1661 143.9934 -223.1430
-78.3448 186.0490 -223.1430 128.0907
If you do no get this result, first compute the matrix elements
of \(P_0(z) = 1\) in the coupled basis: you should get
the identity matrix, because the basis is orthonormal.
Kinetic energy matrix elements
For the expression see exercise 1b of week 5. For the same
basis as above, it is (in cm\(^{-1}\))
T =
0.00000 0.00000 0.00000 0.00000
0.00000 4.02637 0.00000 0.00000
0.00000 0.00000 12.07912 0.00000
0.00000 0.00000 0.00000 24.15823
Tabulate your results in cm\(^{-1}\) as \(E_i-V_\mathrm{min}\).
With the small basis given above, I find
i, E(i) [in cm-1]
0 16.731
1 69.838
2 172.756
3 836.254
Study the convergence of the results
Since this calculation is variational, the energies can not go up in a larger basis.
Report
Instructions for writing you report are given on the main page.
The potential is a one-dimensional cut of the Ar-CO potential from this paper:
J. Chem. Phys. 121, 4691 (2004).
Since I made a one-dimensional cut through the potential, the bending frequency will be higher than in the paper.
HTML5-validator
Last updated: 12-Mar-2021, by Gerrit C. Groenenboom.