Report 2 - Joshua Cooper and Federico Coin
!pip install numba==0.57.0
import numpy as np
from numpy import math
import matplotlib.pyplot as plt
import scipy.special as sp
from numba import njit
from scipy.integrate import quad
from scipy import linalg
from sympy import binomial
import matplotlib.colors as colors
from matplotlib.ticker import MaxNLocator
import warnings
warnings.filterwarnings('ignore')
from IPython.display import Markdown as md
# code for printing matrices
def mprint(matrix):
for row in matrix:
print(' '.join(['{:< 12.2}'.format(x) for x in row]))
return
Question 1
Question 2
input_1
m_p = 938.2720 # masses in Mev
m_n = 939.65653
m = m_p*m_n/(m_p+m_n) # reduced mass
# kinetic energy constant (hbar^2)/(2*m) (MeV*fm)
hbarc = 197.3269804
h_m = hbarc**2 / (2*m)
pi = 4*np.arctan(1)
# B coefficient of the gaussian integral of hermite polynomials
def B_n(n):
values = np.empty(n+1)
values[0] = np.sqrt(pi)
if n>0:
for i in range(1,n+1):
values[i] = values[i-1]*2*i
return values[n]
# Normalization
def Normalisation(n, nu):
value = 1 / np.sqrt(np.sqrt(pi/(32*nu**3)) * 2**(2*n+1) * math.factorial(2*n+1) )
return value
# u_n function
def u_n(r, n, nu):
n_sqrt = np.sqrt(2*nu)
value = Normalisation(n,nu)/n_sqrt * np.exp(-nu*r*r)*(HermitePolynomial(n_sqrt*r,2*n+1))
return value
# Hermite Polynomials
def HermitePolynomial(x, n):
values = np.empty(n+1)
values[0] = 1
values[1] = 2*x
if n>1:
for i in range(2,n+1):
values[i] = 2*x*values[i-1]-2*(i-1)*values[i-2]
return values[n]
ck1 = quad(lambda x: u_n(x, 7, 0.3)**2, 0, np.inf)
print(f'The result of the numerical integration is: {ck1[0]:.1f}.')
# Create x values and prepare y
x = np.linspace(0.01, 6, num=int(6/0.1))
y = np.zeros_like(x)
# Create figure and adjust its properties
fig, ax = plt.subplots(figsize=(10.5, 5.5))
ax.set_title('Basis functions for nu = 0.3')
plt.style.use('ggplot') # Use a different style for the plot
# Calculate y values and plot for each 'i'
for i, linestyle in zip(range(4), ['-', '--', '-.', ':']):
y = [u_n(x_val, i, 0.3)/x_val for x_val in x]
ax.plot(x, y, label=f'n = {i}', linestyle=linestyle) # use different line styles
# Add zero lines and legend
ax.axhline(y=0, color='black', linewidth=.5)
ax.axvline(x=0, color='black', linewidth=.5)
ax.legend(loc='best')
# Set x-label
ax.set_xlabel('r [fm]')
# Show grid
ax.grid(True)
# Change figure background color
fig.patch.set_facecolor('xkcd:light grey')
# Show the plot
plt.show()
# calculating the second derivative
def second_derivative_u_n(x, n, nu):
delta_x = 1e-6 / nu
return (u_n(x + delta_x, n, nu) - 2 * u_n(x, n, nu) + u_n(x - delta_x, n, nu)) / delta_x**2
# defining kinetic energy matrix
def calculate_kinetic_energy(m, n, nu):
nsq, value = np.sqrt(2*nu), Normalisation(n,nu)*Normalisation(m,nu)/np.sqrt(2*nu)*h_m
return value * (
B_n(2 * n + 1) / 4 + (2 * m + 1)*(2 * n + 1)*B_n(2 * n) if n == m else
-n*(2 * n + 1) * B_n(2 * n - 1) if n == (m + 1) else
-m*(2 * m + 1) * B_n(2 * n + 1) if n == (m - 1) else
0
)
# function to compare analytical and numerical results
def checking_kinetic_energy(m, n, nu):
lib_val = quad(lambda x: u_n(x,m, nu)*(-h_m)*second_derivative_u_n(x, n, nu),0,+np.inf)[0]
print(f'For Ekin_{m}_{n}\nOur result: {calculate_kinetic_energy(m,n,nu):.3f}\nNumerical result: {lib_val:.3f}')
checking_kinetic_energy(3, 3, 0.9)
# potential energy constants (MeV)
constants = {'v_r': 200, 'v_t': -178, 'v_s': -91.85}
# width of the potentials (fm^(-2))
widths = {'k_r': 1.487, 'k_t': 0.639, 'k_s': 0.465}
def calculate_temp(i, m, n, g):
return 2**(2*m+1-2*i) * g**(2*(2*(m-i)+1)) * (g**2-1)**(2*i+n-m) * math.exp(
math.log(math.factorial(2*m+1)) + math.log(math.factorial(2*n+1)) - math.log(math.factorial(i)) -\
math.log(math.factorial(2*(m-i)+1)) - math.log(math.factorial(n-m+i)))
def calculate_sum(m, n, g):
i_min = max(0, m - n)
return sum(calculate_temp(i, m, n, g) for i in range(i_min, m+1))
# potential energy matrix elements.
def potential_energy_elements(n ,m, v, k, nu):
g = 1 / np.sqrt(1 + k/(2*nu))
potential_factor = v/(4*nu) * Normalisation(n, nu) * Normalisation(m, nu) * np.sqrt(math.pi / (2*nu + k))
sum_result = calculate_sum(m, n, g)
return potential_factor * sum_result
def verify_potential_calculation(m, n, v, k, nu):
lib_val = quad(lambda x: u_n(x, m, nu) * v * np.exp(-k * x * x) * u_n(x, n, nu), 0, +np.inf)[0]
our_val = potential_energy_elements(m, n, v, k, nu)
print(f'We consider VR_{m}_{n}\npotential_energy_elements yields: {potential_energy_elements(m, n, v, k, nu):.3f}\nThe numerical result is: {lib_val:.3f}')
verify_potential_calculation(2, 5, constants['v_r'], widths['k_r'], 3)
The code below constructs the Hamiltonian matrix from the kinetic and potential energies.
# Defining the size of the Hamiltonian matrix
M = 4
# Allocate space for the Hamiltonian matrix
H = np.zeros((M, M))
# Define the nu parameter
nu = 0.3
# Fill the Hamiltonian matrix
for i in range(M):
for j in range(i, M): # The Hamiltonian is symmetric, so we only need to calculate elements above the diagonal
kinetic = calculate_kinetic_energy(i, j, nu)
potential = potential_energy_elements(i, j, constants['v_r'], widths['k_r'], nu) # using central potential part
H[i, j] = kinetic + potential
if i != j:
H[j, i] = H[i, j] # Use symmetry to fill in the lower triangle
print(H)
Question 3
def compute_ground_state_energy(N_m, nu):
def fill_matrices(V, potential, width):
for i in range(N_m):
for j in range(i+1):
value = potential_energy_elements(i, j, potential, width, nu)
V[i, j] = value
V[j, i] = value
# Create and fill potential and kinetic energy matrices
V_r_1 = np.zeros((N_m, N_m))
fill_matrices(V_r_1, constants['v_r'], widths['k_r'])
V_t_1 = np.zeros((N_m, N_m))
fill_matrices(V_t_1, constants['v_t'], widths['k_t'])
V_s_1 = np.zeros((N_m, N_m))
fill_matrices(V_s_1, constants['v_s'], widths['k_s'])
K_1 = np.zeros((N_m, N_m))
for i in range(N_m):
for j in range(i+1):
value = calculate_kinetic_energy(i, j, nu)
K_1[i, j] = value
K_1[j, i] = value
# Define and calculate eigenvalues for singlet and triplet blocks
singletHamiltonian = V_r_1 + V_s_1 + K_1
singlet_eigenvalues, singlet_eigenvectors = linalg.eigh(singletHamiltonian)
triplet = V_r_1 + V_t_1 + K_1
triplet_eigenvalues, triplet_eigenvectors = linalg.eigh(triplet)
# Find ground state energy and associated eigenvector
ground_state_energy = min(singlet_eigenvalues[0], triplet_eigenvalues[0])
if ground_state_energy == singlet_eigenvalues[0]:
return [ground_state_energy, singlet_eigenvectors[:, 0]]
else:
return [ground_state_energy, triplet_eigenvectors[:, 0]]
num_basis_functions = range(1,40+1)
groundStateEnergies = []
for n in num_basis_functions:
groundStateEnergies.append(compute_ground_state_energy(n,1/8)[0])
from scipy.optimize import curve_fit
# Exponential function
def expFunc(x, a, b, c):
return a * np.exp(-b * x) + c
def plot_ground_state_energies(num_basis_functions, groundStateEnergies):
if len(num_basis_functions) != len(groundStateEnergies):
raise ValueError(f"num_basis_functions and groundStateEnergies must have the same number of elements, but have {len(num_basis_functions)} and {len(groundStateEnergies)} elements respectively")
fig, ax = plt.subplots(figsize=(10.5, 5.5))
ax.plot(num_basis_functions, groundStateEnergies, 'rx')
ax.grid(True)
ax.set_xlabel('$N_{max}$')
ax.set_ylabel('Ground state energy [MeV]')
ax.set_title('Convergence check as a function of Nmax')
# Fitting an exponential function
popt, pcov = curve_fit(expFunc, num_basis_functions, groundStateEnergies, p0=(1, 1e-6, 1))
ax.plot(num_basis_functions, expFunc(num_basis_functions, *popt), 'b--')
ax.text(18, -0.4, 'groundStateEnergies = {:.4} MeV for N = {}'.format(groundStateEnergies[-1], num_basis_functions[-1]),
size=16,
bbox=dict(boxstyle="round", ec=(1., 0.5, 0.5), fc=(1., 0.8, 0.8)))
plot_ground_state_energies(num_basis_functions, groundStateEnergies)
As N increases, we observe a trend towards convergence. Essentially, this implies that the accuracy of our results improves with larger values of N.
The ground state energy derived from this process is negative, characterized by a total spin of one and a total isospin of zero.
These results align with our initial assumptions: the Hamiltonian does not affect isospin, which always maintains a singlet (antisymmetric) state. Given the symmetry of the spatial component of the wavefunction ( 𝑙=0), the spin is required to exist in a triplet (symmetric) state, ensuring the overall state is antisymmetric.
At this point, we're in a position to graphically represent the radial solution of our wavefunction.
def calculate_density(eigenvector, N, nu, x_values):
density = np.zeros(len(x_values))
for n, eigenv in enumerate(eigenvector[:N], start=1):
density += [eigenv * (u_n(x, n - 1, nu) / x) for x in x_values]
return density ** 2
# Get the ground state eigenvector
eigenvector_ground_state = compute_ground_state_energy(40, 1/8)[1]
# Define the x range
distance_array = np.linspace(0.01, 8, num=800)
# Compute the density
density = calculate_density(eigenvector_ground_state, 40, 1/8, distance_array)
# Plot
fig, ax = plt.subplots(figsize=(10.5, 5.5))
ax.plot(distance_array, density)
ax.grid(True)
ax.set_title('Density')
ax.set_ylabel('$|\Psi(r)|^2$')
ax.set_xlabel('r [fm]')
plt.show()
The goal here is to calculate the expected radius for the Ground State. Initially, we need to compute the matrix ⟨𝑛|𝑟^2|𝑚⟩. An in-depth analytical evaluation gives us:
⟨𝑛|𝑟^2|𝑚⟩ is obtained as 1/(𝜈^(√(𝜋) 2^(𝑛+𝑚+2) √((2𝑛+1)!(2𝑚+1)!)) [𝛿_(𝑛,𝑚)(1/2 𝐵_(2𝑛+1)+2(2𝑛+1)(2𝑚+1)𝐵_(2𝑛))+𝛿_(𝑛,𝑚+1)+2𝑛(2𝑛+1)𝐵_(2𝑛−1)+𝛿_(𝑛,𝑚-1)−2𝑚(2𝑚+1)𝐵_(2𝑛+1)].
Thus, by taking |𝐺𝑆⟩ to be the Ground State eigenstate, the expected value for the radius, 𝑟_𝑒, can be formulated as:
𝑟_𝑒=1/2 √(⟨𝐺𝑆|𝑟^2|𝐺𝑆⟩)
def compute_factor(i, j, nnn):
return 1 / (nnn * np.sqrt(np.pi) * 2 ** (i + j + 2) *
np.sqrt(float(math.factorial(2 * i + 1)) * float(math.factorial(2 * j + 1))))
def compute_radius_squared_matrix_element(i, j, factor):
if i == j:
return factor * (0.5 * B_n(2 * i + 1) + 2 * (2 * i + 1) * (2 * j + 1) * B_n(2 * i))
elif i == (j + 1):
return factor * 2 * i * (2 * i + 1) * B_n(2 * i - 1)
elif i == (j - 1):
return factor * 2 * j * (2 * j + 1) * B_n(2 * i + 1)
else:
return 0
def construct_radius_squared_matrix(N_m, nnn):
radius_squared_matrix = np.zeros(shape=(N_m, N_m))
for i in range(N_m):
for j in range(max(0, i - 1), min(N_m - 1, i + 1) + 1):
factor = compute_factor(i, j, nnn)
radius_squared_matrix[i, j] = compute_radius_squared_matrix_element(i, j, factor)
return radius_squared_matrix
eigenVector = compute_ground_state_energy(40, 1/8)[1]
radius_squared_matrix = construct_radius_squared_matrix(40, 1/8)
radius_squared_expectation = eigenVector @ radius_squared_matrix @ eigenVector
radius_expectation = np.sqrt(radius_squared_expectation)
print("The expectation value of the radius of the GS wavefunction is: {:.3} fm".format(radius_expectation/2))
Question 4
We can now investigate convergence across various values of 𝜈. If the width of the Gaussians significantly diverges from the true size of the deuteron, we anticipate convergence will only be achieved for larger 𝑁 values. This phenomenon arises due to the increased complexity in representing a function of specific width using components that are misaligned in size. As evidenced in the visual representation below, convergence is achieved only within the regions depicted in dark blue. This distinctively illustrates that the closer the Gaussian's width approximates the bound state's size, the lower the requisite value for 𝑁 to achieve convergence.
N_values = np.arange(10, 40, 1)
nn_values = np.arange(0.01, 0.5, 0.01)
values = np.zeros((nn_values.size, N_values.size))
for idx_N, N in np.ndenumerate(N_values):
for idx_nn, nn in np.ndenumerate(nn_values):
values[idx_nn, idx_N] = compute_ground_state_energy(N, nn)[0]
from mpl_toolkits.mplot3d import Axes3D
N_values_mesh, nn_values_mesh = np.meshgrid(N_values, nn_values)
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(N_values_mesh, nn_values_mesh, values, cmap='viridis')
fig.colorbar(surf, shrink=0.5, aspect=5)
ax.set_title('groundStateEnergies varying N and nu')
ax.set_xlabel('N')
ax.set_ylabel('nu')
ax.set_zlabel('groundStateEnergies')
plt.show()
def calc_ground_state_energy(matrix_size, nu):
vr_matrix = np.zeros((matrix_size, matrix_size))
vs_matrix = np.zeros((matrix_size, matrix_size))
kinetic_energy_matrix = np.zeros((matrix_size, matrix_size))
for i in range(matrix_size):
for j in range(i + 1):
vr_matrix[i, j] = vr_matrix[j, i] = potential_energy_elements(i, j, constants['v_r'], widths['k_r'], nu)
vs_matrix[i, j] = vs_matrix[j, i] = potential_energy_elements(i, j, constants['v_s'], widths['k_s'], nu)
kinetic_energy_matrix[i, j] = kinetic_energy_matrix[j, i] = calculate_kinetic_energy(i, j, nu)
total_matrix = vr_matrix + vs_matrix + kinetic_energy_matrix
eigenvalues, eigenvectors = linalg.eigh(total_matrix)
return eigenvalues[0], eigenvectors[:, 0]
N_values = range(1, 41)
ground_state_energies = []
for N in N_values:
energy, _ = calc_ground_state_energy(N, 1/100)
ground_state_energies.append(energy)
eigenvect_nn = calc_ground_state_energy(40, 1/100)[1]
print(f'Ground state energy of a n-n system, for nu=0.01, N=40, is: {ground_state_energies[-1]:.3f} MeV')