Commit fa11d9db authored by Ivo Filot's avatar Ivo Filot
Browse files

Tweaking 2D system

parent 7cfd5706
import numpy as np
import matplotlib.pyplot as plt
import scipy
def calculate_mos_1d(N, orbital='s'):
# construct Hamiltonian matrix
......@@ -71,6 +72,8 @@ def construct_hamiltonian_matrix_1d_p(N):
def construct_hamiltonian_matrix_2d(N=32, orbital='s'):
if orbital is 's':
return construct_hamiltonian_matrix_2d_s(N)
if orbital is 'p':
return construct_hamiltonian_matrix_2d_pxyz(N)
if orbital is 'dxy':
return construct_hamiltonian_matrix_2d_dxy(N)
return None
......@@ -92,6 +95,11 @@ def construct_hamiltonian_matrix_2d_s(N):
return m
def construct_hamiltonian_matrix_2d_pxyz(N):
return -1.0 * scipy.linalg.block_diag(construct_hamiltonian_matrix_2d_s(N),
construct_hamiltonian_matrix_2d_s(N),
construct_hamiltonian_matrix_2d_s(N))
def construct_hamiltonian_matrix_2d_dxy(N):
m = np.zeros((N*N,N*N))
for i in range(N):
......@@ -102,20 +110,10 @@ def construct_hamiltonian_matrix_2d_dxy(N):
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
......@@ -159,7 +157,8 @@ def plot_coefficients_1d(coefficients, idx):
plt.ylabel('$c_{i}$')
plt.show()
def plot_coefficients_2d(coefficients, idx):
def plot_coefficients_2d(coefficients, idx, contour=True,
interpolation='quadric'):
coeff = coefficients.transpose()[idx]
n = np.sqrt(len(coeff))
......@@ -168,11 +167,17 @@ def plot_coefficients_2d(coefficients, idx):
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)
maxval = 2.0 / n
if contour:
contours = plt.contour(x, y, z, colors='black', zorder=2,
levels=np.linspace(-maxval, maxval, 13))
plt.clabel(contours, inline=True, fontsize=8)
plt.imshow(z, extent=[0, n, 0, n], origin='lower',
interpolation='quadric', vmin=-maxval, vmax=maxval)
else:
plt.imshow(z, extent=[0, n, 0, n], origin='lower',
interpolation='None', vmin=-maxval, vmax=maxval)
ax.set_aspect(1)
plt.xlabel('$i$')
plt.ylabel('$j$')
......
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='p')
# 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
......@@ -11,4 +11,7 @@ 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
huckel.plot_coefficients_2d(v, i)
# plot result for highest lying orbital
huckel.plot_coefficients_2d(v, -1, contour=False, interpolation='None')
\ 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