Python file selfconsistent

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()