Skip to content

Solutions for bonds and spectra

Warm-up exercises

1. $$ \begin{pmatrix} E_0 & -t & -\tilde{t} \ -t & E_0 & -t \ -\tilde{t} & -t & E_0 \end{pmatrix} $$

2. $$ \begin{pmatrix} E_0 & -t & 0 & 0 & 0 & 0\ -t & E_0 & -t & 0 & 0 & 0\ 0 & -t & E_0 & -t & 0 & 0\ 0 & 0 & -t & E_0 & -t & 0\ 0 & 0 & 0 & -t & E_0 & -t\ 0 & 0 & 0 & 0 & -t & E_0 \end{pmatrix} $$

3.Look at the sign of the force \(F = -\frac{d U(r)}{dr}\).

Exercise 1: linear triatomic molecule

  1. In 1D, there are two normal modes and in 3D there are 4 normal modes
  2. $$ \begin{cases} m \ddot{u}_1 & = -\kappa(u_1-u_2)\ M \ddot{u}_2 & = -\kappa(2u_2-u_1-u_3)\ m \ddot{u}_3 & = -\kappa(u_3-u_2) \end{cases} $$ Where \(m\) is the mass of the oxygen atoms and \(M\) the mass of the carbon atom.
  3. \[ \omega = \sqrt{\frac{\kappa}{m}} \]
  4. \[ \frac{u_1}{u_2} = -\frac{M}{2m} \]
  5. \[ \omega = \sqrt{\frac{\kappa(2m+M)}{mM}} \]

Exercise 2: Lennard-Jones potential

  1. See lecture+internet
  2. The equilibrium position is \(r_0 = 2^{1/6}\sigma\). The energy at the inter atomic distance \(r_0\) is given by: $$ U(r_0) = -\epsilon $$
  3. $$ U(r) = -\epsilon + \frac{\kappa}{2}(r-r_0)^2 $$ Where \(\kappa = \frac{72\epsilon}{2^{1/3}\sigma^2}\)
  4. The ground state energy is given by $$ E_0 = -\epsilon+\frac{1}{2}\hbar\sqrt{\frac{2\kappa}{m}} $$ And the breaking energy is given by $$ E_{break} = \epsilon - \frac{1}{2}\hbar\sqrt{\frac{2\kappa}{m}} $$
  5. Distance from which \(U(r)\) becomes anharmonic: $$ r_{anharmonic} = \frac{6}{7}2^{1/6}\sigma $$ Number of phonons that fit in the potential before it becomes anharmonic $$ n = \frac{1}{49\hbar}\sqrt{\frac{\kappa m}{2}}-\frac{1}{2} $$

Exercise 3: Numerical simulation

1.

2.

import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt

# Creating function that initializes the spring matrix
def initial_mat(N):
    K = np.zeros((N,N), dtype = int)
    b = np.ones(N-1)
    np.fill_diagonal(K, 2); np.fill_diagonal(K[1:], -b); np.fill_diagonal(K[:, 1:], -b) 
    K[0,0] = K[-1,-1] = 1
    omega = np.sqrt(np.abs(LA.eigvalsh(K)))

    # outputting a histogram of the data
    plt.figure()
    plt.hist(omega, bins = 30)
    plt.xlabel("$\omega$")
    plt.ylabel("Number of levels per eigenfrequency")
    plt.title('Number of levels per eigenfrequencies for N = %d' %N)

# Running the code
initial_mat(5)
initial_mat(200)

png

png

3.

import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt

# Creating mass matrix
def mass_matrix(N, m1, m2):
    M = np.zeros((N,N), dtype = int)
    j = np.linspace(0, N-1, N)
    for j in range(N):
        if j%2 ==0:
            M[j, j] = m1
        else:
            M[j, j] = m2
    return M

# Creating function that initializes the spring matrix
def initial_mat(N,M):
    K = np.zeros((N,N), dtype = int)
    b = np.ones(N-1)
    np.fill_diagonal(K, 2); np.fill_diagonal(K[1:], -b); np.fill_diagonal(K[:, 1:], -b) 
    K[0,0] = K[-1,-1] = 1
    MK = np.dot(LA.inv(M), K)
    omega = np.sqrt(np.abs(LA.eigvalsh(MK)))

    # outputting a histogram of the data
    plt.figure()
    plt.hist(omega, bins = 30)
    plt.xlabel("$\omega$")
    plt.ylabel("Number of levels per eigenfrequency")

# Defining variables
N = 200
m1 = 1; m2 = 4

# Running the code
M = mass_matrix(N, m1, m2)
initial_mat(N, M)

png

Where in this simulation, every even numbered atom in the chain has 4 times the mass of every odd numbered atom

4.

# Defining constants
kappa_1 = 1
kappa_2 = 2

# Creating function that initializes the spring matrix
def initial_mat(N, kappa_1, kappa_2):
    K = np.zeros((N,N), dtype = int)
    b = np.ones(N-1)
    c = np.zeros(N-1)
    idx = np.linspace(0, N-2, N-1)

    # Create pattern of zero's and ones
    b[idx % 2 == 0] = 0
    c[idx % 2 == 0] = 1

    # Diagnonal
    np.fill_diagonal(K, kappa_1+kappa_2)

    # Off diagnonal
    np.fill_diagonal(K[1:], -c*kappa_1-b*kappa_2)
    np.fill_diagonal(K[:, 1:], -c*kappa_1-b*kappa_2)  


    # Setting up initial and last values
    K[0,0] = kappa_1
    if (N % 2) == 0:
        K[-1,-1] = kappa_1
    else:
        K[-1,-1] = kappa_2

    # Calculating the dispersion relation
    omega = np.sqrt(np.abs(LA.eigvalsh(K)))

    # outputting a histogram of the data
    plt.figure()
    plt.hist(omega, bins = 30)
    plt.xlabel("$\omega$")
    plt.ylabel("Number of levels per eigenfrequency")
    plt.title('Number of levels per eigenfrequencies for N = %d' %N)


# Running the code
initial_mat(200, kappa_1, kappa_2)

png