Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- from scipy import integrate
- from scipy.linalg import eigh_tridiagonal
- import numpy as np
- import matplotlib.pyplot as pyplot
- from scipy.sparse import diags
- from scipy.linalg import eigh_tridiagonal
- import timeit
- import time
- start_time = time.time()
- # list of physical constants
- const_1 = 0.0380998 # h^2/(2*m) in nm^2 eV
- const_2 = 1.43996 # e^2/(4*pi*e0) in nm eV where e0 is the permittivity of free space
- r0 = 0.0529177 # bohr radius in nm
- h = 6.62606896e-34 # planks constant
- c = 299792458 # speed of light in ms^-1
- hc = 1239.8419 # planck constant * speed of light in eV nm
- alpha = 0.01 # change to 0.01 for task 1
- # parameters for size of matrices and maximum r
- n = 12000
- r_max = 2
- # this is the analytic solution to the potential integral from r to infinity as a function of r for r > 0
- analytic_sol = lambda r: (const_2)*(((const_2)/(2*const_1))**alpha) * (1/(alpha-1)) * (r**(alpha - 1))
- # lambda function for the electrostatic force as a function of r
- force = lambda r: (-const_2)*(1/(r**2))*(( const_2*r/(2*const_1))**alpha)
- # create empty list to add to later
- r_list = [] # list of r positions
- v_list = [] # list of potentials corresponding to r_list
- a_list = [] # list of analytic solutions of potentials corresponding to r list
- diff_list = [] # absolute values of the differences between analytic and numerical solutions
- # function will be used to create a list of positions r from
- # also returns the distance between each of these positions, del_r
- def create_r_list(n, end):
- lst = []
- increment = end/n
- start = end - n*increment
- for i in range(n):
- lst.append(start + i*increment)
- lst.append(end)
- del lst[0]
- return lst , increment
- # function returns the lowest two eigenvalues of a matrix
- def find_evals(mat):
- e_vals, e_vecs = eigh_tridiagonal(np.diag(mat), np.diag(mat, k=1))
- e0 = min(e_vals)
- delete = np.array([e0])
- e1 = min(np.setdiff1d(e_vals, delete))
- return e0, e1
- # function uses scipy to integrate electrostatic force from r to infinity to find electrostatic potential
- def integrate_force(f,r):
- return integrate.quad(f, r, np.inf)[0]
- # create list of r ranging from 0.01nm to 1 nm
- r_list, del_r = create_r_list(100, 1)
- # adds to v_list the calculated potentials from the r in r list, and a_list the corresponding analytic solution
- # adds to diff_list the absolute difference between corresponding values in these two lists
- for r in r_list:
- v = integrate_force(force,r)
- a = analytic_sol(r)
- v_list.append(v)
- a_list.append(a)
- diff_list.append(abs(v - a))
- # plotting potential vs r graph
- pyplot.title('Potential V(r) vs r')
- pyplot.xlabel('r,nm')
- pyplot.ylabel('V(r), eV')
- pyplot.plot(r_list, v_list)
- pyplot.show()
- # prints the biggest relative difference, diff, between the exact potential and calculated potential
- print("The least accurate potential is off by " + str(max(diff_list)) + " eV")
- alpha_list = [0.0, 0.01]
- for a in alpha_list:
- alpha = a
- # now resetting r and v lists as we are using larger n and r_max to calculate accurate energy levels
- r_list = []
- v_list = []
- r_list, del_r = create_r_list(n, r_max)
- # adds to v list the calculated potentials from the r in r list, and a_list the corresponding analytic solution
- for r in r_list:
- v = integrate_force(force,r)
- v_list.append(v)
- del_mat = (1/(del_r)**2)*diags([1, -2, 1], [-1, 0, 1], shape=(n , n )).toarray()
- v_mat = np.diag(v_list)
- h_mat = -const_1*del_mat + v_mat
- e0, e1 = find_evals(h_mat)
- print("For alpha = " + str(alpha) + ":")
- print("The first energy level is: " + str(e0))
- print("While the second energy level is: " +str(e1))
- # creates list of n+1 numbers from start number to end number
- def create_alpha_list(n, start, end):
- lst = []
- increment = (end - start)/n
- for i in range(n):
- lst.append(start + i*increment)
- lst.append(end)
- return lst
- # create list of values of alpha from 0 to 0.01
- alpha_list = create_alpha_list(15, 0, 0.01)
- # empty list to add energy difference between first and second energy levels, calculated from corresponding alpha in alpha_list
- diff_list = []
- for a in alpha_list:
- alpha = a
- v_list = []
- for r in r_list:
- v = integrate_force(force,r)
- v_list.append(v)
- v_mat = np.diag(v_list)
- h_mat = -const_1*del_mat + v_mat
- e0, e1 = find_evals(h_mat)
- diff_list.append(abs(e0-e1))
- now = time.time()
- print(start_time - now)
- # plotting alpha vs energy level difference
- pyplot.title('Alpha vs Energy Difference Between First and Second Energy Level')
- pyplot.xlabel('Alpha')
- pyplot.ylabel('Energy Level Differencs, eV')
- pyplot.plot(alpha_list, diff_list)
- pyplot.show()
- # calculating the energy difference we desire between level 1 and 2
- wavelength = 121.503
- energy_aim = hc/wavelength
- # calculating the lower and upper bounds of acceptable energy differences between level 1 and 2
- wavelength_lower = 121.503 - 0.001
- energy_aim_upper = hc/wavelength_lower
- wavelength_upper = 121.503 + 0.001
- energy_aim_lower = hc/wavelength_upper
- if diff_list[0] < energy_aim:
- print("pheeeeeewwwwwiiieeeee")
- else:
- print("uh oh")
- # subtracting the
- diff_list_new = [d - energy_aim for d in diff_list]
- energy_aim_upper_new = energy_aim_upper - energy_aim_upper
- energy_aim_lower_new = energy_aim_lower - energy_aim_upper
- energy_aim_new = energy_aim - energy_aim_upper
- def secant(x1,x2,fx1,fx2):
- x3 = x2 - fx2*( (x2 - x1) / (fx2 - fx1))
- return x3
- # finding a sensible upper bound to bracket the root
- b = 0.0
- diff = 0.0
- for d in diff_list_new:
- if d < energy_aim_upper_new:
- diff = d
- index = diff_list_new.index(d)
- b = alpha_list[index + 1]
- diff = diff_list_new[index + 1]
- check = 0
- count = 0
- x1 = alpha_list[0]
- x2 = b
- fx1 = diff_list_new[0]
- fx2 = diff
- while check == 0:
- print(x1,x2,fx1,fx2)
- new_x1 = x2
- new_fx1 = fx2
- x2 = secant(x1, x2, fx1, fx2)
- print("x2", x2)
- alpha = x2
- v_list = []
- for r in r_list:
- v = integrate_force(force,r)
- v_list.append(v)
- v_mat = np.diag(v_list)
- h_mat = -const_1*del_mat + v_mat
- e_vals, e_vecs = eigh_tridiagonal(np.diag(h_mat), np.diag(h_mat, k=1))
- e0, e1 = find_evals(h_mat)
- fx2 = abs(e0 - e1) - energy_aim_upper
- count += 1
- print(count)
- print(x1)
- print(x2)
- if abs(x2 - x1) < 0.000001 and fx2 < energy_aim_upper_new:
- check = 1
- x1 = new_x1
- fx1 = new_fx1
- print("yooooooooo")
- now = time.time()
- print(start_time - now)
- print(" ")
- print(x2)
- print(fx2)
- print(fx2 + energy_aim_upper, energy_aim_lower, energy_aim, energy_aim_upper)
- print(fx2 +energy_aim_upper < energy_aim_upper)
- alpha = x2
- v_list = []
- for r in r_list:
- v = integrate_force(force,r)
- v_list.append(v)
- v_mat = np.diag(v_list)
- h_mat = -const_1*del_mat + v_mat
- e_vals, e_vecs = eigh_tridiagonal(np.diag(h_mat), np.diag(h_mat, k=1))
- e0, e1 = find_evals(h_mat)
- fx2 = abs(e0 - e1)
- print(fx2, energy_aim_upper)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement