In [1]:

    

%pylab inline
import sys, os
import QuickPlot
import matplotlib.pylab as pylab
pylab.rcParams['figure.figsize'] = 8, 6  # that's default image size for this interactive session
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
from ipywidgets.widgets import interactive, fixed, interact
%config InlineBackend.figure_format = 'retina'


    

    




Populating the interactive namespace from numpy and matplotlib






    




/Users/Zack/anaconda/envs/conda36/lib/python3.6/site-packages/matplotlib/__init__.py:1405: UserWarning: 
This call to matplotlib.use() has no effect because the backend has already
been chosen; matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
or matplotlib.backends is imported for the first time.

  warnings.warn(_use_error_msg)






    

In [10]:

    

hr = genfromtxt('NiO_hr.dat', skip_header=22)
print(hr.shape)


    

    




(17856, 7)

In [11]:

    

# This reads the hamiltonian elements for a specific unit cell coordinate.  
# Wannier90 reports them by unit cell (x,y,z, i,j) where (1,0,0, 1,1) is the overlap of
# wavefunction 1 in the 0,0,0 cell with wavefunction 1 translated by one unit cell in the x direction.
def GetOneCell(x,y,z):
    hrtrim = hr[(hr[:,0]==x) & (hr[:,1]==y) & (hr[:,2]==z)]
    hrtrim = hrtrim[:,3:]
    H = np.zeros((8,8), dtype='complex')
    for r in hrtrim[:]:
        i = int(r[0])-1
        j = int(r[1])-1
        overlap = complex(r[2], r[3])
        H[i, j] = overlap
    return H


    

In [12]:

    

# Print the eigenvalues for just the consideration of one unit cell.
H = GetOneCell(0,0,0)
print(eig(H)[0])


    

    




[ 10.82645182 -5.75837398e-16j  10.82778016 -4.28947938e-16j
  10.82732174 -8.46326285e-16j  13.82136642 +4.25342641e-17j
  13.82057484 +2.94782951e-16j  13.57947358 +1.38536426e-15j
  13.58115549 -1.63569424e-15j  13.58040296 -1.00492796e-15j]

In [13]:

    

# Now do the same but considering the 6 nearest neighbors.
H = GetOneCell(0,0,0)

H += GetOneCell(1,0,0)
H += GetOneCell(0,1,0)
H += GetOneCell(0,0,1)
H += GetOneCell(-1,0,0)
H += GetOneCell(0,-1,0)
H += GetOneCell(0,0,-1)

print(eig(H)[0])


    

    




[ 10.27226110 +9.96714859e-16j  10.27084844 +7.00582207e-16j
  11.29528462 +8.60729761e-16j  14.49484033 -8.61192853e-17j
  14.49292364 -5.26749628e-16j  13.26279285 -1.45946758e-15j
  13.26655856 -3.65687579e-17j  13.83063147 +1.32766894e-15j]

In [14]:

    

# Now build it using every element.
H = np.zeros((8,8), dtype='complex')
for r in hr[:]:
    x = int(r[0])
    y = int(r[1])
    z = int(r[2])
    i = int(r[3])-1
    j = int(r[4])-1
    overlap = complex(r[5], r[6])
    H[i, j] += overlap

print(eig(H)[0])


    

    




[ 13.32768716 +1.23611329e-16j  13.33084012 +2.79505415e-16j
  13.03345671 -8.98853128e-16j  12.96499531 +3.52657773e-16j
  12.96381012 +8.11242317e-16j  11.50521747 -3.80358771e-16j
  11.50449025 -3.68450604e-16j  11.50476785 +1.92047443e-16j]

In [47]:

    

# Print out the energy levels from greatest to least, 
# along with <psi_m|psi_n> so we can see what orbitals
# they represent.
E = linalg.eig(H)[0]
psi = linalg.eig(H)[1]

for i in range(len(E)):
    print(f'Energy = {np.abs(E[i]):0.2f} eV')
    print(f'psi=(', end='')
    normval=np.abs(psi[i]).sum()
    for j in range(len(E)):
        print(f'{np.abs(psi[i][j])/normval:0.2f}, ', end='')
    print(')')


    

    




Energy = 13.33 eV
psi=(0.60, 0.39, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, )
Energy = 13.33 eV
psi=(0.00, 0.00, 0.10, 0.45, 0.44, 0.00, 0.00, 0.00, )
Energy = 13.03 eV
psi=(0.00, 0.00, 0.58, 0.14, 0.27, 0.00, 0.00, 0.00, )
Energy = 12.96 eV
psi=(0.39, 0.60, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, )
Energy = 12.96 eV
psi=(0.00, 0.00, 0.25, 0.39, 0.35, 0.00, 0.00, 0.00, )
Energy = 11.51 eV
psi=(0.00, 0.00, 0.00, 0.00, 0.00, 0.30, 0.40, 0.30, )
Energy = 11.50 eV
psi=(0.00, 0.00, 0.00, 0.00, 0.00, 0.28, 0.38, 0.34, )
Energy = 11.50 eV
psi=(0.00, 0.00, 0.00, 0.00, 0.00, 0.43, 0.19, 0.38, )

In [ ]: