Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import matplotlib.pyplot as plt
- import numpy as np
- from tqdm import tqdm
- # Constants
- m = 0.5109989461e6
- h = 0.19732697
- C = 1.44e-3
- Z = 2
- N = 100
- L = 0.5e-3
- dr = (L-0.1e-3)/N
- R = np.linspace(0.1e-3,L,N)
- Uscf = np.zeros(N) # Inital Guess
- fig, (ax1, ax2) = plt.subplots(1, 2,figsize=(10,5))
- for blah in range(25):
- constants = -(h**2)/(2*m)
- laplacian = -2 * np.eye(N) + np.eye(N, k=1) + np.eye(N, k=-1)
- laplacian /= dr ** 2
- columb = -Z * C * (1/R)
- H = constants * laplacian + np.diag(columb) + np.diag(Uscf)
- vals, vecs = np.linalg.eig(H)
- indicies = np.argsort(vals)
- vals = vals[indicies]
- vecs = vecs[indicies]
- f = vecs[:,0]
- sigma = np.abs(f)**2
- for i,r in enumerate(R):
- Uscf[i] = C * (np.sum(sigma[:i])/r + np.sum(sigma[i+1:]))
- ax1.plot(R, sigma/R)
- ax1.set_title('E = ' + str(vals[0]))
- ax1.set_ylabel('Psi^Psi [nm^{-1}]')
- ax1.set_xlabel('distance [nm]')
- ax2.plot(R, columb, label='Columb')
- ax2.plot(R, Uscf, linestyle='--', label='SCF')
- ax2.set_ylabel('energy [eV]')
- ax2.set_xlabel('distance [nm]')
- ax2.legend()
- fig.tight_layout()
- plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement