Commit 7cfd5706 authored by Ivo Filot's avatar Ivo Filot
Browse files

Adding 2D structures

parent 76d47dd1
# Byte-compiled / optimized / DLL files
__pycache__/
*.py[cod]
*$py.class
# C extensions
*.so
# Distribution / packaging
.Python
build/
develop-eggs/
dist/
downloads/
eggs/
.eggs/
lib/
lib64/
parts/
sdist/
var/
wheels/
share/python-wheels/
*.egg-info/
.installed.cfg
*.egg
MANIFEST
# PyInstaller
# Usually these files are written by a python script from a template
# before PyInstaller builds the exe, so as to inject date/other infos into it.
*.manifest
*.spec
# Installer logs
pip-log.txt
pip-delete-this-directory.txt
# Unit test / coverage reports
htmlcov/
.tox/
.nox/
.coverage
.coverage.*
.cache
nosetests.xml
coverage.xml
*.cover
.hypothesis/
.pytest_cache/
# Translations
*.mo
*.pot
# Django stuff:
*.log
local_settings.py
db.sqlite3
# Flask stuff:
instance/
.webassets-cache
# Scrapy stuff:
.scrapy
# Sphinx documentation
docs/_build/
# PyBuilder
target/
# Jupyter Notebook
.ipynb_checkpoints
# IPython
profile_default/
ipython_config.py
# pyenv
.python-version
# celery beat schedule file
celerybeat-schedule
# SageMath parsed files
*.sage.py
# Environments
.env
.venv
env/
venv/
ENV/
env.bak/
venv.bak/
# Spyder project settings
.spyderproject
.spyproject
# Rope project settings
.ropeproject
# mkdocs documentation
/site
# mypy
.mypy_cache/
.dmypy.json
dmypy.json
# Pyre type checker
.pyre/
# cache
__pycache__
*.pyc
\ No newline at end of file
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
def calculate_mos_1d(N, orbital='s'):
# construct Hamiltonian matrix
m = construct_hamiltonian_matrix(N, orbital)
return m
# 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 calculate_mos(N):
def calculate_mos_2d(N, orbital='s'):
# construct Hamiltonian matrix
m = construct_hamiltonian_matrix(N)
m = construct_hamiltonian_matrix_2d(N, orbital)
# obtain eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eigh(m)
......@@ -27,6 +29,100 @@ def calculate_mos(N):
return -eigenvalues, eigenvectors
##############################################################################
# Hamiltonian matrices for 1D systems
##############################################################################
def construct_hamiltonian_matrix(N=1000, orbital='s'):
if orbital is 's':
return construct_hamiltonian_matrix_1d_s(N)
if orbital is 'p':
return construct_hamiltonian_matrix_1d_p(N)
return None
def construct_hamiltonian_matrix_1d_s(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 construct_hamiltonian_matrix_1d_p(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
##############################################################################
# Hamiltonian matrices for 2D systems
##############################################################################
def construct_hamiltonian_matrix_2d(N=32, orbital='s'):
if orbital is 's':
return construct_hamiltonian_matrix_2d_s(N)
if orbital is 'dxy':
return construct_hamiltonian_matrix_2d_dxy(N)
return None
def construct_hamiltonian_matrix_2d_s(N):
m = np.zeros((N*N,N*N))
for i in range(N):
for j in range(N):
atid0 = i * N + j
atid1 = np.mod(i+1, N) * N + j
atid2 = np.mod(i-1, N) * N + j
atid3 = i * N + np.mod(j+1, N)
atid4 = i * N + np.mod(j-1, N)
m[atid0, atid1] = 1
m[atid0, atid2] = 1
m[atid0, atid3] = 1
m[atid0, atid4] = 1
return m
def construct_hamiltonian_matrix_2d_dxy(N):
m = np.zeros((N*N,N*N))
for i in range(N):
for j in range(N):
atid0 = i * N + j
atid1 = np.mod(i+1, N) * N + j
atid2 = np.mod(i-1, N) * N + j
atid3 = i * N + np.mod(j+1, N)
atid4 = i * N + np.mod(j-1, N)
atid5 = np.mod(i+1, N) * N + np.mod(j+1, N)
atid6 = np.mod(i-1, N) * N + np.mod(j+1, N)
atid7 = np.mod(i+1, N) * N + np.mod(j-1, N)
atid8 = np.mod(i-1, N) * N + np.mod(j-1, N)
m[atid0, atid1] = -1
m[atid0, atid2] = -1
m[atid0, atid3] = -1
m[atid0, atid4] = -1
m[atid0, atid5] = 1
m[atid0, atid6] = 1
m[atid0, atid7] = 1
m[atid0, atid8] = 1
return m
##############################################################################
# Plot functions
##############################################################################
def plot_histogram(e, perc):
nn = int(len(e)*perc)
avg = np.average(e[0:nn])
......@@ -37,21 +133,24 @@ def plot_histogram(e, perc):
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.hist([set1, set2], edgecolor='black', color=['grey', '0.9'],
bins=bins, zorder=2, label=['Occupied', '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.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])
def plot_coefficients_1d(coefficients, idx):
coeff = coefficients.transpose()[idx]
Norb = len(coeff)
plt.figure()
xrange = np.arange(0, len(coefficients.transpose()[idx]))
plt.plot(xrange, coefficients.transpose()[idx], 'o')
xrange = np.arange(0, len(coeff))
plt.plot(xrange, coeff, 'o')
maxcoeff = np.sqrt(2.0 / (Norb + 1.0)) * 1.1
plt.ylim(-maxcoeff, maxcoeff)
......@@ -59,6 +158,29 @@ def plot_coefficients(coefficients, idx):
plt.xlabel('$i$')
plt.ylabel('$c_{i}$')
plt.show()
def plot_coefficients_2d(coefficients, idx):
coeff = coefficients.transpose()[idx]
n = np.sqrt(len(coeff))
x = np.arange(0, int(n))
y = np.arange(0, int(n))
z = np.reshape(np.array(coeff), (int(n), int(n)))
fig, ax = plt.subplots()
contours = plt.contour(x, y, z, colors='black', zorder=2,
levels=np.linspace(-0.06, 0.06, 13))
plt.clabel(contours, inline=True, fontsize=8)
plt.imshow(z, extent=[0, n, 0, n], origin='lower',
interpolation='quadric', vmin=-0.06, vmax=0.06)
ax.set_aspect(1)
plt.xlabel('$i$')
plt.ylabel('$j$')
plt.colorbar();
plt.xticks(np.linspace(0, n, 5))
plt.yticks(np.linspace(0, n, 5))
plt.gcf().subplots_adjust(bottom=0.15)
plt.show()
def plot_mo(eigenvalues):
plt.figure()
......
import calc.phuckel as huckel
# set number of orbitals
Norb = 1000
# calculate energies and eigenvectors
e, v = huckel.calculate_mos_1d(Norb, orbital='p')
# plot a histogram of the band energies
huckel.plot_histogram(e, 0.50)
# plot the coefficients of the lowest and highest lying orbital
huckel.plot_coefficients_1d(v, 0)
huckel.plot_coefficients_1d(v, -1)
\ No newline at end of file
import calc.phuckel as huckel
# set number of orbitals
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
# calculate energies and eigenvectors
e, v = huckel.calculate_mos_1d(Norb, orbital='s')
# plot a histogram of the band energies
huckel.plot_histogram(e, 0.50)
# plot the coefficients of the four lowest lying orbitals
for i in range(0,4):
huckel.plot_coefficients_1d(v, i)
\ No newline at end of file
import calc.phuckel as huckel
# set number of orbitals
Norb = 32 # will be 32 x 32 array
# calculate energies and eigenvectors
e, v = huckel.calculate_mos_2d(Norb, orbital='s')
# plot a histogram of the band energies
huckel.plot_histogram(e, 0.50)
# plot four degenerate states z=1-4
for i in range(1,5):
huckel.plot_coefficients_2d(v, i)
\ 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