Advertisement
Guest User

Untitled

a guest
Oct 13th, 2019
173
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 7.21 KB | None | 0 0
  1. from scipy import integrate
  2. from scipy.linalg import eigh_tridiagonal
  3. import numpy as np
  4. import matplotlib.pyplot as pyplot
  5. from scipy.sparse import diags
  6. from scipy.linalg import eigh_tridiagonal
  7. import timeit
  8. import time
  9.  
  10. start_time = time.time()
  11.  
  12. # list of physical constants
  13. const_1 = 0.0380998 # h^2/(2*m) in nm^2 eV
  14. const_2 = 1.43996 # e^2/(4*pi*e0) in nm eV where e0 is the permittivity of free space
  15. r0 = 0.0529177 # bohr radius in nm
  16. h = 6.62606896e-34 # planks constant
  17. c = 299792458 # speed of light in ms^-1
  18. hc = 1239.8419 # planck constant * speed of light in eV nm
  19. alpha = 0.01 # change to 0.01 for task 1
  20.  
  21. # parameters for size of matrices and maximum r
  22. n = 12000
  23. r_max = 2
  24.  
  25. # this is the analytic solution to the potential integral from r to infinity as a function of r for r > 0
  26. analytic_sol = lambda r: (const_2)*(((const_2)/(2*const_1))**alpha) * (1/(alpha-1)) * (r**(alpha - 1))
  27.  
  28. # lambda function for the electrostatic force as a function of r
  29. force = lambda r: (-const_2)*(1/(r**2))*(( const_2*r/(2*const_1))**alpha)
  30.  
  31. # create empty list to add to later
  32. r_list = [] # list of r positions
  33. v_list = [] # list of potentials corresponding to r_list
  34. a_list = [] # list of analytic solutions of potentials corresponding to r list
  35. diff_list = [] # absolute values of the differences between analytic and numerical solutions
  36.  
  37. # function will be used to create a list of positions r from
  38. # also returns the distance between each of these positions, del_r
  39. def create_r_list(n, end):
  40. lst = []
  41. increment = end/n
  42. start = end - n*increment
  43. for i in range(n):
  44. lst.append(start + i*increment)
  45. lst.append(end)
  46. del lst[0]
  47. return lst , increment
  48.  
  49. # function returns the lowest two eigenvalues of a matrix
  50. def find_evals(mat):
  51. e_vals, e_vecs = eigh_tridiagonal(np.diag(mat), np.diag(mat, k=1))
  52. e0 = min(e_vals)
  53. delete = np.array([e0])
  54. e1 = min(np.setdiff1d(e_vals, delete))
  55.  
  56. return e0, e1
  57.  
  58.  
  59. # function uses scipy to integrate electrostatic force from r to infinity to find electrostatic potential
  60. def integrate_force(f,r):
  61. return integrate.quad(f, r, np.inf)[0]
  62.  
  63. # create list of r ranging from 0.01nm to 1 nm
  64. r_list, del_r = create_r_list(100, 1)
  65.  
  66. # adds to v_list the calculated potentials from the r in r list, and a_list the corresponding analytic solution
  67. # adds to diff_list the absolute difference between corresponding values in these two lists
  68. for r in r_list:
  69. v = integrate_force(force,r)
  70. a = analytic_sol(r)
  71. v_list.append(v)
  72. a_list.append(a)
  73. diff_list.append(abs(v - a))
  74.  
  75. # plotting potential vs r graph
  76. pyplot.title('Potential V(r) vs r')
  77. pyplot.xlabel('r,nm')
  78. pyplot.ylabel('V(r), eV')
  79. pyplot.plot(r_list, v_list)
  80. pyplot.show()
  81.  
  82. # prints the biggest relative difference, diff, between the exact potential and calculated potential
  83. print("The least accurate potential is off by " + str(max(diff_list)) + " eV")
  84.  
  85. alpha_list = [0.0, 0.01]
  86.  
  87. for a in alpha_list:
  88. alpha = a
  89. # now resetting r and v lists as we are using larger n and r_max to calculate accurate energy levels
  90. r_list = []
  91. v_list = []
  92. r_list, del_r = create_r_list(n, r_max)
  93.  
  94. # adds to v list the calculated potentials from the r in r list, and a_list the corresponding analytic solution
  95. for r in r_list:
  96. v = integrate_force(force,r)
  97. v_list.append(v)
  98.  
  99. del_mat = (1/(del_r)**2)*diags([1, -2, 1], [-1, 0, 1], shape=(n , n )).toarray()
  100. v_mat = np.diag(v_list)
  101. h_mat = -const_1*del_mat + v_mat
  102.  
  103. e0, e1 = find_evals(h_mat)
  104. print("For alpha = " + str(alpha) + ":")
  105. print("The first energy level is: " + str(e0))
  106. print("While the second energy level is: " +str(e1))
  107.  
  108.  
  109.  
  110.  
  111. # creates list of n+1 numbers from start number to end number
  112. def create_alpha_list(n, start, end):
  113. lst = []
  114. increment = (end - start)/n
  115. for i in range(n):
  116. lst.append(start + i*increment)
  117. lst.append(end)
  118. return lst
  119.  
  120. # create list of values of alpha from 0 to 0.01
  121. alpha_list = create_alpha_list(15, 0, 0.01)
  122. # empty list to add energy difference between first and second energy levels, calculated from corresponding alpha in alpha_list
  123. diff_list = []
  124.  
  125.  
  126. for a in alpha_list:
  127. alpha = a
  128. v_list = []
  129. for r in r_list:
  130. v = integrate_force(force,r)
  131. v_list.append(v)
  132.  
  133. v_mat = np.diag(v_list)
  134. h_mat = -const_1*del_mat + v_mat
  135.  
  136. e0, e1 = find_evals(h_mat)
  137. diff_list.append(abs(e0-e1))
  138. now = time.time()
  139. print(start_time - now)
  140.  
  141.  
  142. # plotting alpha vs energy level difference
  143. pyplot.title('Alpha vs Energy Difference Between First and Second Energy Level')
  144. pyplot.xlabel('Alpha')
  145. pyplot.ylabel('Energy Level Differencs, eV')
  146. pyplot.plot(alpha_list, diff_list)
  147. pyplot.show()
  148.  
  149.  
  150. # calculating the energy difference we desire between level 1 and 2
  151. wavelength = 121.503
  152. energy_aim = hc/wavelength
  153.  
  154. # calculating the lower and upper bounds of acceptable energy differences between level 1 and 2
  155. wavelength_lower = 121.503 - 0.001
  156. energy_aim_upper = hc/wavelength_lower
  157. wavelength_upper = 121.503 + 0.001
  158. energy_aim_lower = hc/wavelength_upper
  159.  
  160.  
  161. if diff_list[0] < energy_aim:
  162. print("pheeeeeewwwwwiiieeeee")
  163. else:
  164. print("uh oh")
  165.  
  166. # subtracting the
  167. diff_list_new = [d - energy_aim for d in diff_list]
  168. energy_aim_upper_new = energy_aim_upper - energy_aim_upper
  169. energy_aim_lower_new = energy_aim_lower - energy_aim_upper
  170. energy_aim_new = energy_aim - energy_aim_upper
  171.  
  172. def secant(x1,x2,fx1,fx2):
  173. x3 = x2 - fx2*( (x2 - x1) / (fx2 - fx1))
  174. return x3
  175.  
  176. # finding a sensible upper bound to bracket the root
  177. b = 0.0
  178. diff = 0.0
  179. for d in diff_list_new:
  180. if d < energy_aim_upper_new:
  181. diff = d
  182. index = diff_list_new.index(d)
  183. b = alpha_list[index + 1]
  184. diff = diff_list_new[index + 1]
  185.  
  186.  
  187. check = 0
  188. count = 0
  189.  
  190. x1 = alpha_list[0]
  191. x2 = b
  192. fx1 = diff_list_new[0]
  193. fx2 = diff
  194.  
  195.  
  196. while check == 0:
  197. print(x1,x2,fx1,fx2)
  198. new_x1 = x2
  199. new_fx1 = fx2
  200. x2 = secant(x1, x2, fx1, fx2)
  201. print("x2", x2)
  202.  
  203. alpha = x2
  204. v_list = []
  205. for r in r_list:
  206. v = integrate_force(force,r)
  207. v_list.append(v)
  208. v_mat = np.diag(v_list)
  209. h_mat = -const_1*del_mat + v_mat
  210. e_vals, e_vecs = eigh_tridiagonal(np.diag(h_mat), np.diag(h_mat, k=1))
  211. e0, e1 = find_evals(h_mat)
  212.  
  213. fx2 = abs(e0 - e1) - energy_aim_upper
  214. count += 1
  215. print(count)
  216. print(x1)
  217. print(x2)
  218. if abs(x2 - x1) < 0.000001 and fx2 < energy_aim_upper_new:
  219. check = 1
  220. x1 = new_x1
  221. fx1 = new_fx1
  222. print("yooooooooo")
  223. now = time.time()
  224. print(start_time - now)
  225.  
  226. print(" ")
  227. print(x2)
  228. print(fx2)
  229. print(fx2 + energy_aim_upper, energy_aim_lower, energy_aim, energy_aim_upper)
  230. print(fx2 +energy_aim_upper < energy_aim_upper)
  231.  
  232. alpha = x2
  233. v_list = []
  234. for r in r_list:
  235. v = integrate_force(force,r)
  236. v_list.append(v)
  237. v_mat = np.diag(v_list)
  238. h_mat = -const_1*del_mat + v_mat
  239. e_vals, e_vecs = eigh_tridiagonal(np.diag(h_mat), np.diag(h_mat, k=1))
  240. e0, e1 = find_evals(h_mat)
  241.  
  242. fx2 = abs(e0 - e1)
  243. print(fx2, energy_aim_upper)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement