COMPLETE CODE

Homework 1

Students: Federico Coin, Joshua Cooper

!pip install numba==0.56.4

import numpy as np import matplotlib.pyplot as plt from numba import njit import math import scipy.special as sp from scipy.signal import find_peaks

Step 1

# Define Lennard-Johnes potential def pot(x, s): #depends on the variables r (x) and sigma (s) return 4 * eps * ((s / x) ** 12 - (s / x) ** 6)

from scipy.constants import c,hbar m_hbar2 = 939 * 1e9 / (0.197 * 1e3 * 1e4) ** 2 eps = 5.99 #epsilon parameter pi = 4*np.arctan(1)

The following block contains the Numerov's method.

# Define K^2 def K2(E, x, l, s): return 2 * m_hbar2 * (E - pot(x, s)) - l * (l + 1) / x ** 2 # Numerov's method def numerical_solution(E, l, x, s): y = np.empty(len(x)) h = abs(x[1] - x[0]) y[0] = np.exp(-(s / 2 / x[0]) ** (5)) y[1] = np.exp(-(s / 2 / x[1]) ** (5)) for i in range(np.size(x) - 2): y[i + 2] = 1 / (1 + 1 / 12 * h ** 2 * K2(E, x[i + 2], l, s)) * (y[i + 1] * (2 - 5 / 6 * h ** 2 * K2(E, x[i + 1], l, s)) - y[i] * (1 + 1 / 12 * h ** 2 * K2(E, x[i], l, s))) return y

E_0 = 1.5 sigma_0 = 2.0

# Mesh definition x_min = 0.5 # starting point of the mesh x_max = 30 # maximum value of x step = 0.01 # step value l_max = 7 # angular momentum x_0 = np.arange(x_min, x_max, step)

# Curves plot fig = plt.figure() fig.set_size_inches(10.5, 5.5) plt.title("Solution of the SE for different values of $l$") for l in range(0, l_max): y = numerical_solution(E_0, l, x_0, sigma_0) plt.plot(x_0, y / np.amax(y[int((15 - x_min) / step):]), label="l = " + str(l)) plt.xlabel("r [Å]") plt.ylabel("Wave functions") plt.grid("on") plt.legend() plt.show()

Step 2

def bessel_spherical_jl(l, x): j = np.zeros(len(x)) for i, xx in enumerate(x): j[i] = np.sqrt(np.pi / (2 * xx)) * sp.jv(l + 0.5, xx) return j def bessel_spherical_nl(l, x): n = np.zeros(len(x)) for i, xx in enumerate(x): n[i] = np.sqrt(np.pi / (2 * xx)) * sp.yv(l + 0.5, xx) return n x_0 = np.linspace(0, 20, 200) ll = np.array([1, 2]) jl_data = [bessel_spherical_jl(l, x_0) for l in ll] nl_data = [bessel_spherical_nl(l, x_0) for l in ll] fig, axs = plt.subplots(2, 2, figsize=(12, 8)) fig.suptitle('Spherical Bessel and Neumann Functions') styles = ['-', '--', ':'] for i in range(2): for j in range(2): if j == 0: axs[i, j].plot(x_0, jl_data[i], styles[i], label='Custom', linewidth=1) axs[i, j].plot(x_0, sp.spherical_jn(ll[i], x_0), '*r', label='Scipy', linewidth=0.1) axs[i, j].set_title('Spherical Bessel, l={}'.format(ll[i])) else: axs[i, j].plot(x_0, nl_data[i], styles[i], label='Custom', linewidth=1) axs[i, j].plot(x_0, sp.spherical_yn(ll[i], x_0), '*r', label='Scipy', linewidth=0.1) axs[i, j].set_title('Spherical Neumann, l={}'.format(ll[i])) axs[i, j].set_ylim([-2, 1]) axs[i, j].legend(loc='best') plt.show()

The graphs above clearly show that the results are consistent with the implementations from the scipy library.

def bessel(lmax,r): v1 = np.cos(r)/r v2 = np.sin(r)/r j = np.array([v1,v2]) n = np.array([v2,-v1]) l0 = 0 while l0<lmax: j = np.append(j,((2*l0+1)/r*j[-1]-j[-2])) n = np.append(n,((2*l0+1)/r*n[-1]-n[-2])) l0 += 1 return [j,n] def phase_shift_tan(energy, l, r1, r2, x_vals, y_vals): k = np.sqrt(2 * energy * m_hbar2) idx_r1, idx_r2 = int((r1 - x_min) / step), int((r2 - x_min) / step) x1, x2 = x_vals[idx_r1], x_vals[idx_r2] beta = (y_vals[idx_r1] * x2) / (y_vals[idx_r2] * x1) j1, n1 = bessel(l, k * x1) j2, n2 = bessel(l, k * x2) return (beta * j2[-1] - j1[-1]) / (beta * n2[-1] - n1[-1]) x_vals = np.linspace(0.5, 100.01, num=int((100.01-0.5)/0.01)) energy, l_value = 1.5, 5 y_vals = numerical_solution(energy, l_value, x_vals, 2.5) tan_deltas = [phase_shift_tan(energy, l_value, r, r + 1, x_vals, y_vals) for r in np.arange(5, 61, 0.01)] #calculate_shift_tangent_v2 fig, ax = plt.subplots() ax.plot(np.arange(5, 61, 0.01), tan_deltas, linestyle='--', marker='o', markersize=4, color='red') ax.set_xlabel('$r_1$', fontsize=14) ax.set_ylabel('$tan(\delta_l)$', fontsize=14) ax.set_title(f'Convergence of the phase shift', fontsize=16) ax.text(20, 0.08, 'We kept the distance between $r_1,r_2$ fixed to 1', fontsize=12) ax.grid(True) plt.show()

Step 3

Here we define some functions that we will be using for step 3. The 'scattering_coefficient' function takes the energy, phase shift and angular momenta as inputs to calculate the partial cross section term, while the other functions calculate the conversion between energy and velocity. These functions are used in order to calculate the cross section.

@njit def velocity_to_energy(velocity_value): return 0.5 * (velocity_value / c) ** 2 * 939 * 1e9 @njit def energy_to_velocity(energy_value): return np.sqrt(2 * energy_value / (939 * 1e9)) * c @njit def scattering_coefficient(energy, l_val, shift_tan): return 4 * np.pi / (2 * energy * m_hbar2) * (2 * l_val + 1) * (shift_tan ** 2) / (1 + shift_tan ** 2)

The block below contains a function that computes the values of the total cross section computed on the values of energy of the array of energies values. The inputs to the function is the and the maximum value of angular momentum l, value of sigma, and an array of energy values.

def compute_cross_section(sigma_val, energy_list, r1_val, r2_val, lmax): cross_section_matrix = np.empty((lmax + 1, len(energy_list))) for idx, energy in enumerate(energy_list): for l in range(lmax + 1): y_vals = numerical_solution(energy, l, x, sigma_val) cross_section_matrix[l, idx] = scattering_coefficient(energy, l, phase_shift_tan(energy, l, r1_val, r2_val, x, y_vals)) return np.sum(cross_section_matrix, axis=0)

# Defining the Energy Mesh E_min = 0.2 E_max = 3.5 E_step = 0.0075 energies = np.arange(E_min,E_max + E_step,E_step) # Angular momentum l_max = 8 # Position mesh x = np.arange(x_min,x_max,step) # Position value for the phaseshift r_1 = 25 r_2 = 26

Next we plot as shown below.

figure, axis = plt.subplots(constrained_layout=True) figure.set_size_inches(10.5, 5.5) # Set plot style plt.style.use('ggplot') axis.plot(energy_to_velocity(energies), compute_cross_section(sigma_0, energies, r_1, r_2, lmax=6), label='l_max = 6', linestyle='--', marker='o', markersize=4) axis.plot(energy_to_velocity(energies), compute_cross_section(sigma_0, energies, r_1, r_2, lmax=8), label='l_max = 8', linestyle='-.', marker='s', markersize=4) axis.plot(energy_to_velocity(energies), compute_cross_section(sigma_0, energies, r_1, r_2, lmax=4), label='l_max = 4', linestyle=':', marker='D', markersize=4) #axis.plot(energy_to_velocity(energies), compute_cross_section(sigma_0, energies, r_1, r_2, lmax=10), label='l_max = 10', linestyle=':', marker='D', markersize=4) axis.grid(True) secondary_axis = axis.secondary_xaxis('top', functions=(velocity_to_energy, energy_to_velocity)) secondary_axis.set_xlabel("E [meV]") axis.set_xlabel('Relative speed [m/s]') axis.set_ylabel('Cross section [Å^2]') axis.set_title('H-Kr cross section, σ = {:.3f} Å'.format(sigma_0)) # axis.plot(exp_data1[:, 0], exp_data1[:, 1], 'r.', label='Experimental data, Grabbed') plt.semilogx() plt.legend() plt.show()

Step 4

In this section we used the experimental values reported in Table I of Toennis' paper, which are the positions of the total cross section peaks in energy.

(DISCLAIMER: the following block can take up to two hours to run)

def calculate_ratios_diffs_and_chi_squared(sigma_try, e_try, experimental_energies, experimental_error, r_1, r_2, l_max): calculated_cross_section = compute_cross_section(sigma_try, e_try, r_1, r_2, lmax=l_max) idxs = find_peaks(calculated_cross_section)[0] energy_peaks = e_try[idxs][:3] chi_squared = np.sum((energy_peaks - experimental_energies) ** 2 / experimental_error ** 2) difference = energy_peaks - experimental_energies ratio = difference / experimental_error return ratio, difference, chi_squared # Experimental values of the positions of the total cross section peaks in energy experimental_energies = np.array([0.50, 1.59, 2.94]) experimental_error = np.array([0.02, 0.06, 0.12]) e_try = np.arange(E_min, E_max, 0.0075) # After trying out several values, we used this interval sigma_vals = np.arange(3.07, 3.13, 0.0005) results = [calculate_ratios_diffs_and_chi_squared(sigma_try, e_try, experimental_energies, experimental_error, r_1, r_2, l_max) for sigma_try in sigma_vals] ratio, difference, chi_squared = zip(*results) optimal_sigma = sigma_vals[np.argmin(chi_squared)] print('σ = {:.4f} +- {} Å^2'.format(optimal_sigma, 0.0005))

cut_off_value = 1.5 store = np.where((abs(np.array(ratio)[:, 0]) <= cut_off_value) & (abs(np.array(ratio)[:, 1]) <= cut_off_value) & (abs(np.array(ratio)[:, 2]) <= cut_off_value)) uncertainty_sigma = np.amax(sigma_vals[store] - optimal_sigma) print('σ = {:0.3f} +- {:0.3f} Å^2'.format(optimal_sigma, uncertainty_sigma))

Below, the code plots the results.

#import matplotlib.pyplot as plt #from utilities import compute_cross_section, e2v, v2e cross_section_final = compute_cross_section(optimal_sigma, energies, r_1, r_2, lmax=l_max) fig, ax = plt.subplots(constrained_layout=True) fig.set_size_inches(10.5, 5.5) plt.style.use('seaborn-darkgrid') ax.set_facecolor('#F0F0F0') secondary_axis = ax.secondary_xaxis('top', functions=(energy_to_velocity, velocity_to_energy)) ax.set_title('H-Kr cross section, $σ = {:.3f} Å$'.format(optimal_sigma), fontsize=18, fontweight='bold') secondary_axis.set_xlabel('Relative speed $[m/s]$', fontsize=14) ax.set_xlabel("E [meV]", fontsize=14) ax.set_ylabel('Cross section $[Å^2]$', fontsize=14) ax.plot(energies, cross_section_final, label='Computed values', color='black', linestyle='-', linewidth=1) exp_energy_data = [{'energy': experimental_energies[i], 'color': color, 'error': experimental_error[i]} for i, color in enumerate(['blue', 'orange', 'cyan'])] for data in exp_energy_data: ax.vlines(data['energy'], 100, 600, linestyles=':', color=data['color'], label=f' $E = {data["energy"]} meV$', linewidth=2) plt.axvspan(data['energy'] - data['error'], data['energy'] + data['error'], alpha=0.25, color=data['color']) plt.legend(fontsize=12) plt.show()

As demonstrated in this plot, the 3 peaks that we have calculated agree with the the values given in the paper which is further confirmed with the uncertainty regions.