6_bonds_and_spectra_solutions.md

Solutions for bonds and spectra

Warm-up exercises

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

$$\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

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)

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

# 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)