• API
• FAQ
• Tools
• Archive
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
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):
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 < 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
191. x2 = b
192. fx1 = diff_list_new
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.

Top