Commit 76d47dd1 authored by Ivo Filot's avatar Ivo Filot
Browse files

Adding script for 1d periodic system

parent 83c4c216
import numpy as np
import matplotlib.pyplot as plt
def construct_hamiltonian_matrix(N):
m = np.zeros((N,N))
for i in range(N):
if i < N-1:
m[i,i+1] = 1
if i>0:
m[i,i-1] = 1
m[0,N-1] = 1
m[N-1,0] = 1
return m
def calculate_mos(N):
# construct Hamiltonian matrix
m = construct_hamiltonian_matrix(N)
# obtain eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eigh(m)
# sort eigenvalues and eigenvectors
idx = eigenvalues.argsort()[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:,idx]
return -eigenvalues, eigenvectors
def plot_histogram(e, perc):
nn = int(len(e)*perc)
avg = np.average(e[0:nn])
var = np.var(e[0:nn])
set1 = e[0:nn]
set2 =e[nn:len(e)]
bins=np.histogram(np.hstack((set1,set2)), bins=16)[1] #get the bin edges
plt.figure()
plt.hist([set1, set2], edgecolor='black', color=['grey', '0.9'], bins=bins, zorder=2, label=['Occupied', 'Unoccupied'], stacked=True)
#plt.hist(set2, edgecolor='black', color='0.9', bins=bins, zorder=2, label='Unoccupied', stacked=True)
plt.grid(linestyle='--', color='grey', zorder=1)
plt.legend()
plt.xlabel('Energy $E$')
plt.ylabel('Occurence')
plt.title('$\\mu = %4.2f \\beta, \\sigma = %4.2f \\beta$' % (avg, np.sqrt(var)))
plt.gcf().subplots_adjust(left=0.20, bottom=0.20)
plt.show()
def plot_coefficients(coefficients, idx):
Norb = len(coefficients.transpose()[idx])
plt.figure()
xrange = np.arange(0, len(coefficients.transpose()[idx]))
plt.plot(xrange, coefficients.transpose()[idx], 'o')
maxcoeff = np.sqrt(2.0 / (Norb + 1.0)) * 1.1
plt.ylim(-maxcoeff, maxcoeff)
plt.grid(linestyle='--', color='0.7')
plt.xlabel('$i$')
plt.ylabel('$c_{i}$')
plt.show()
def plot_mo(eigenvalues):
plt.figure()
plt.grid(linestyle='--', color='0.7')
plt.ylim((-2.5,2.5))
plt.plot(np.zeros(len(eigenvalues)), eigenvalues, 'o', alpha=0.5)
plt.show()
\ No newline at end of file
import calc.phuckel as huckel
Norb = 1000
e, v = huckel.calculate_mos(Norb)
huckel.plot_histogram(e, 0.25)
huckel.plot_mo(e)
huckel.plot_coefficients(v, 1)
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment