Entry 1947

# Polycrystalline Mobility Edge Model (Street, Northrup, Salleo)
# Jeremy Smith
# Imperial College London

from numpy import *
from scipy import array
from scipy.optimize import leastsq
import os

# Constants
e = 1.60217646e-19
k = 1.3806503e-23
A = 1.06e56

# Function definitions
def listMean(x):
	n, mean = len(x), 0
	for a in x:
		mean = mean + log(a)
	mean = mean / float(n)
	return mean

def num_integrate(y, x, a, b):
	"""Numerical integration using trapezium rule between x[a] and x[b]"""
	trap = [0]
	for i in range(len(y)-1):
		t = 0.5*(y[i] + y[i+1])*(x[i+1] - x[i])
		trap.append(t)
	intg = 0
	for i in trap[a:b+1]:
		intg += i
	return intg

#Variables for material
Ntot = 1.5e26   #PARAM
eb_meV = 120   #PARAM
eb = eb_meV*e/1000

mu0 = 10.0   #PARAM
mu0B = 1.0   #PARAM

em_meV = 30   #PARAM
em = em_meV*e/1000

#Variables for device
Vgate = [30.0, 40.0, 50.0]
Ci = 2.1e-5
l = 1e-9


#Energy variable
E_eV = [float(x)/1000 for x in range(-500, 1200)]
E = [x*e for x in E_eV]

#Fermi energy starting guess
ef_eV = 0.1
ef = ef_eV*e

#Opening and reading data file

datafile = open(os.path.normpath("C:/Documents and Settings/jnsmith/My Documents/Python scripts/Polyaxial ME model/data-full.txt"),"r")
mu_expt = []
T = []
for line in datafile:
	a,b = line.strip().split("\t")
	T.append(float(a))
	mu_expt.append(float(b))
datafile.close()

#Total sum of squares (to calculate r)
SStot = 0
mu_expt_mean = listMean(mu_expt)
for y in mu_expt:
	SStot += (log(y) - mu_expt_mean)**2

#Fitting to the data
params = [Ntot, eb, em, mu0, mu0B]

def modelResiduals(p, y, tmp, vg):
	err = log(y) - log(mobEdgeModel(p, tmp, vg))
	print "err:", err
	print "params:", p
	print "r:", sqrt(1 - (sum(err**2)/SStot))
	print "\n"
	return err

def mobEdgeModel(p, tmp, vg):
	muEFF = []
	dataLoop = 0
	for v in vg:
		print " Vg =", v
		dataLoop += 1
		#Density of states and mobility models
		beta = (p[0]/(p[1]*A))**2
		dos = [A*sqrt(beta - z) for z in E if z<0] + [(p[0]/abs(p[1]))*exp(-z/abs(p[1])) for z in E if z>=0]
		mobility = [p[3] for z in E if z<0] + [p[4]*exp(-z/abs(p[2])) for z in E if z>=0]
		nTdevice = Ci*v/(e*l)
			
		def fermi(x, x0, t):
			"""Fermi-Dirac distribution"""
			return 1/(1 + exp((x-x0)/(k*t)))

		def residuals(q, y, t, d):
			err = log10(y) - log10(nTcalculator(q, t, d))
			return err
			
		def nTcalculator(q, t, d):
			"""Calculates the total density of charge carriers given the DOS"""
			occupied_dos = []
			for i in range(len(E)):
				occupied_dos.append(d[i]*fermi(-E[i], -q[0], t))
			return num_integrate(occupied_dos, E, 0, len(E))
		
		#Solves the total charge carrier concentration by altering Ef
		tmp_partial = tmp[0:14]    #Selects the correct temperatures for a single vg
		m = []
		for t in tmp_partial:
			ef_fit = leastsq(residuals, [ef], args=(nTdevice, t, dos), maxfev=20)
			#print "   T=", t, "   Ef=", ef_fit[0]/e

			#Mobility calculation
			occupied_dos_mobility = []
			for i in range(len(E)):
				occupied_dos_mobility.append(dos[i]*fermi(-E[i], -ef_fit[0], t)*mobility[i])
			m.append(num_integrate(occupied_dos_mobility, E, 0, len(E))/nTdevice)
		muEFF += m
	return array(muEFF)
	
print "\nTemperature list:"
print T
print "Mobility data:"
print mu_expt
print "Gate voltages:"
print Vgate
print "\n\n"

fitModel = leastsq(modelResiduals, params, args=(mu_expt, T, Vgate), ftol=1.5e-7, xtol=1.5e-7, maxfev=1500, factor=10)
print "\n"
print "Ntot =", fitModel[0][0]
print "Eb =", 1000*fitModel[0][1]/e
print "Em =", 1000*fitModel[0][2]/e
print "mu0 =", fitModel[0][3]
print "mu0B =", fitModel[0][4]
print "Starting parameters: ", [Ntot, eb_meV, em_meV, mu0, mu0B]

print "Complete."
raw_input()