Python file selfconsistent
selfconsistent.py — Python Source, 1Kb
File contents
import sys import numpy as np import pylab as plt from constants import * # Parameters U0 = 0.025 kT = 0.025 mu = 0 ep = 0.2 g1 = 0.005 g2 = 0.005 g = g1 + g2 alphag = 0.9 # C1/C, eq. (1.15) alphad = 0.1 # C2/C, eq. (1.15) # Energy grid NE = 501 E = np.linspace(-1, 1, NE) dE = E[1] - E[0] # Lorentzian Density of states per eV D = g / (2 * np.pi) / (E**2 + (g / 2)**2) print D.sum() * dE D = D / (dE * D.sum()) # normalizing to one print D.sum() * dE # Bias IV = 101 V = np.linspace(0, 6, IV) dV = V[1] - V[0] #dV = 1. / (IV - 1) #V = np.arange(0, + dV, dV) N = np.empty(V.shape) I = np.empty(V.shape) for iV, VV in enumerate(V): Vg = 0 Vd = VV # Vd = 0 # Vg = VV[iV] mu1 = mu mu2 = mu1 - Vd # energies in eV UL = - (alphag * Vg) - (alphad * Vd) U = 0 # self-consistent field dU = 1 while dU > 1e-6: f1 = 1. / (1 + np.exp((E + ep + UL + U - mu1) / kT)) f2 = 1. / (1 + np.exp((E + ep + UL + U - mu2) / kT)) N[iV] = dE * np.sum(D * ((f1 * g1 / g) + (f2 * g2 / g))) Unew = U0 * N[iV] dU = abs(U - Unew) U = U + 0.1 * (Unew - U) ## print "iV, dU=", iV, dU I[iV] = dE * I0 * np.sum(D * (f1 - f2) * g1 * g2 / g) plt.figure() plt.subplot(211) plt.grid() plt.plot(V, I) plt.xlabel('Voltage [V]') plt.ylabel('Current [A]') plt.subplot(212) plt.grid() plt.plot(V, N) plt.xlabel('Voltage [V]') plt.ylabel('# electrons') plt.show()