Computational and Theoretical Chemistry II, Computer assigment 2

CTC2 (NWI-MOL175): Computer assignment 2: variational calculation atom-diatom system

Read the entire assignment before you start coding.

Compute the bending modes of the Ar-CO complex with a variational calculation using a coupled angular momentum basis \( |(jl)JM\rangle \). 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

Matlab and Python functions to compute Clebsch-Gordan coefficients are available here.
This special case of Clebsch-Gordan coefficients may simplify your code: \[ \langle j_1 m_1 j_2 m_2 | 0\; 0 \rangle = \frac{(-1)^{j_1-m_1}}{\sqrt{2j_1+1}} \delta_{j_1,j_2} \delta_{m_1,-m_2}. \]

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: 26-Apr-2022, by Gerrit C. Groenenboom.