

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)

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)

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)