import numpy as np from scipy.special import erf, dawsn from scipy.optimize import fsolve, bisect import scipy sqrt2 = np.sqrt(2) def naplot(e, Delta=None, Lambda=None, eta=None, epsa=None, epsF=None): import pylab as pl if Delta is not None: pl.plot(e, Delta, 'r', lw=2) if Lambda is not None: pl.plot(e, Lambda, 'b', lw=2) if eta is not None: pl.plot(e, eta, 'g', lw=2) if epsa is not None: pl.plot(e, e - epsa, 'k--', lw=2) pl.plot([epsa], [0.0], 'ko') if epsF is not None: pl.axvline(epsF, ls='k--') class NewnsAnderson: def __init__(self, epsk, epsa, vak, n, width=0.03, gk=None): """n is the number of substrate electrons.""" self.epsk = np.array(epsk, float) self.epsa = epsa self.vak = np.array(vak, float) self.n = n self.width = width if gk is None: gk = np.ones(len(epsk), dtype=int) self.gk = gk self.Mg = gk.sum() def plot(self, e, Delta=None, Lambda=None, eta=None, epsa=None, epsF=None): naplot(e, Delta, Lambda, eta, epsa, epsF) def G(self, eps): return 0.5 * (1.0 + erf(eps / (sqrt2 * self.width))) def g(self, eps): A = 1.0 / (np.sqrt(2.0 * np.pi) * self.width) e = np.exp(-0.5 * (eps / self.width)**2) return A * e def gt(self, eps): "evaluates the Hilbert transform of the gaussian" v = sqrt2 / (np.pi * self.width) * dawsn(eps / (sqrt2 * self.width)) return v def get_dos(self, eps): dos = np.zeros_like(eps) for eps_k, g_k in zip(self.epsk, self.gk): dos += self.g(eps - eps_k) * g_k return dos def get_dos_integral(self, eps0): return self.G(eps0 - self.epsk).sum() def get_epsF(self, dn=0.0): if dn == 'homo': return self.epsk.repeat(self.gk)[(self.n - 1) // 2] elif dn == 'lumo': return self.epsk.repeat(self.gk)[self.n // 2] def count(eps): return 2.0 * self.get_dos_integral(eps) - self.n - dn eps0 = (self.epsk[0] + self.epsk[-1]) / 2.0 if self.epsk[-1] - self.epsk[0] < 1e-4: epsF = 0.5 * (self.epsk[0] + self.epsk[-1]) else: epsF = bisect(count, self.epsk[0] - 0.01, self.epsk[-1] + 0.01, xtol=1e-5, #xtol=10.0 * self.width ) #epsF = fsolve(count, epsF, fprime=self.get_dos) return epsF def atan(self, a, b): #+pi - yes? "Returns arctan (a,b)+pi, such that arctan(a/b)=arctan(b,a)" j = scipy.sqrt(-1) v = -j * scipy.log((a + b * j) / (scipy.sqrt(a * a + b * b))) - np.pi return v.real def Delta(self, eps): """The imaginary part of the self-energy""" # pi*sum_k |Vak|**2 g(eps-eps_k) res = np.zeros_like(eps) for eps_k, g, v in zip(self.epsk, self.gk, self.vak): #print type(v), type(g), type(eps), type(eps_k) res += v**2 * g * self.g(eps - eps_k) res *= np.pi return res def Lambda(self, eps): """The real part of the self-energy Sigma""" # sqrt(2)/sigma * sum_k |Vak|**2 F^dav((eps-eps_k)/sqrt(2)*sigma)) res = np.zeros_like(eps) for eps_k, g, v in zip(self.epsk, self.gk, self.vak): res += v**2 * g * self.gt(eps - eps_k) res *= np.pi return res def eta(self, eps, D=None, L=None): """ The induced phase change""" if D is None: D = self.Delta(eps) if L is None: L = self.Lambda(eps) v = self.atan((eps - self.epsa - L), D) return v def build_hamiltonian(self): M = self.gk.sum() + 1 H = np.zeros((M, M)) H[:-1, :-1] = np.diag(self.epsk.repeat(self.gk)) H[-1, -1] = self.epsa vak = self.vak.repeat(self.gk) H[:-1, -1] = vak H[-1, :-1] = vak return H def diagonalize(self): H = self.build_hamiltonian() from scipy.linalg import eigvalsh eigs = eigvalsh(H) n = self.n #print len(eigs), len(eigs[:-1]), len(self.gk) E = eigs.repeat(2)[:n + 1].sum() Ek = self.epsk.repeat(2 * self.gk)[:n].sum() Ea = self.epsa dE = E - Ek - Ea #print n, dE, E, Ek, self.epsk.repeat(2 * self.gk)[n-1] #print E, Ek, Ea return dE def get_adsorption_energy(self, e, eta=None, epsF=None): if eta is None: eta = self.eta(e) de = e[1] - e[0] if epsF is None: epsF = self.get_epsF('homo') etasum = eta[e<=epsF].sum() * de / np.pi dE = 2.0 * etasum + epsF - self.epsa return dE def delta_rho(self, eps, eta=None): "induced changed in density of states : not implemented yet" # d/dx arctan(B[x],x-A[x]) # = # B(x)(1-A'(x))+(A(x)-x)B'(x) #---------------------------- # (x-A(x))^2+B(x)^2 # # however I will do this numerically if eta is None: eta=self.eta(eps) deta=[0.0] for i in range(len(eta)-2): deta_1=(eta[i+2]-eta[i+1])/(eps[i+1]-eps[i]) deta_2=(eta[i+1]-eta[i])/(eps[i+1]-eps[i]) deta.append((deta_1+deta_2)/2.0) deta.append(0) deta=np.array(deta) return -1.0 / np.pi * deta def pdos(self, e, Delta=None, Lambda=None): if Delta is None: Delta = self.Delta(e) if Lambda is None: Lambda = self.Lambda(e) drho = Delta / ((e - self.epsa - Lambda)**2 + Delta**2) / np.pi return drho def get_eta_integral(self, e, eta=None, epsF=None): if eta is None: eta = self.eta(e) if epsF is None: epsF = self.get_epsF('homo') de = e[1] - e[0] integral = eta[e <= epsF].sum() * de / np.pi return integral def main(): x = np.linspace(-10.0, 10.0, 4000) dx = x[1] - x[0] na = NewnsAnderson([-4,-1.94,0.1], -3, [2,1,2], 4, width=0.1 ) print na.get_dos_integral(-1.98) epsF = na.get_epsF() print epsF print na.get_dos(x).sum() * dx print na.get_dos(x[x <= epsF]).sum() * dx emin = na.epsk[0] - 10.0 * na.width emax = epsF + 2.0 egrid = np.linspace(emin, emax, 2048) D = na.Delta(egrid) L = na.Lambda(egrid) eta = na.eta(egrid, D=D, L=L) acc_e = np.cumsum(eta) * (dx / np.pi * 2.0) import pylab as pl P = 3 pl.subplot(P,1,1) pl.plot(egrid, D) pl.plot(egrid, L) pl.subplot(P,1,2) pl.plot(egrid, eta) pl.subplot(P, 1, 3) pl.plot(egrid, acc_e) pl.show() if __name__ == '__main__': main()