SHARE
TWEET

Untitled

a guest Oct 13th, 2019 88 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  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)
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top