molecule.py

File contents

# -*- coding: utf-8 -*-
"""
Created on Thu Jun 25 16:44:46 2015

@author: walt
"""
import numpy as np
import pylab as plt

eps0 = -0.25
eps1 = 0.25
tau = 0.5

eta = 0.25
eta = 0.00001
H = np.array([[eps0, tau], [tau, eps1]])

E = np.linspace(-1, 1, 501)
D = np.empty((len(E), 2))
for i, e in enumerate(E):
    G = np.linalg.inv((e + 1j * eta) * np.eye(2) - H)
    A = 1j * (G - G.conj().T)
    D[i, 0] = A[0, 0].real
    D[i, 1] = A[1, 1].real
D = D.T

plt.plot(E, D[0], label='LDOS(0)', color='r')
plt.plot(E, D[1], label='LDOS(1)', color='b')

# self energy treatment
Sigma = tau**2 / (E - eps1 + 1j * eta)
g = 1.0 / (E + 1j * eta - eps0 - Sigma)
a =  1j * (g - g.conj())
plt.plot(E, a.real, '--', label='g(0)', color='r', lw=2)

plt.legend()
plt.xlabel('energy [eV]')
plt.ylabel('LDOS')
plt.show()