#TODO: BulkGrapheneWithTemperature.Ef_interp auf genauigkeit testen! min/maxEf parametrisieren
import sympy.mpmath
import scipy.optimize
from .common import Constants
import numpy
import math
import scipy.interpolate
class QuantumCapacitanceSelfConsistency:
"""
Solve the quantum capacitance equation self-consistently through
iteration. It works, but the convergence is not thoroughly tested,
thus the class is DEPRECATED. Use QuantumCapacitanceSolver instead.
"""
container=0
charge_operator=0
solver=0
convergence_tol=0
def __init__(self,container,charge_operator,convergence_tol=1e-2,solver=None):
if solver!=None:
self.solver=solver
else:
self.solver,inhomogeneity=container.lu_solver()
self.charge_operator=charge_operator
self.container=container
self.convergence_tol=convergence_tol
def self_consistency_cycle(self,elements):
rectangle_elementnumbers_range=self.container.rectangle_elementnumbers_range()
while True:
inhom=self.container.createinhomogeneity()
solution=self.solver(inhom)
diff=0
for elem in elements:
charge=self.container.charge(solution,self.charge_operator,[elem])
old=elem.potential
elem.potential=elem.fermi_energy-elem.fermi_energy_charge_dependence(charge)/Constants.elem_charge
diff+=math.fabs(elem.potential-old)
print diff
if diff < self.convergence_tol:
break
inhom=self.container.createinhomogeneity()
solution=self.solver(inhom)
return solution
#Solution can also be obtained by looking at the potential value of the elements
[docs]class QuantumCapacitanceSolver:
"""
Calculates the quantum capacitance of a system with a density of states
D(E). D(E,x) should work too (with equation='charge_nonlocal), but hasn't
been tested yet.
Note: QuantumCapacitanceSolver changes the potential property of the
participating elements. The constant property is fermi_energy (i.e.
if you connect an external voltage of 5V, set fermi_energy=5).
Usage:
1) Create object: qcsolver=QuantumCapacitanceSolver(container,solver,
elements,operator)
2) Calculate classical dependency of charge on every possible potential
configuration of elements: qcsolver.refresh_basisvecs()
3) Set potential/fermi_energy as you like (e.g. backgate, sidegates...).
4) Refresh the contribution of the new configuration:
qcsolver.refresh_environment_contrib()
5) Solve:
qcsolution=qcsolver.solve()
6) go to 3) until you are through with all voltages you want to look at.
"""
container=None
solver=None
elements=None
charge_operator=None
__equation=None
__nonlocal_charge_potential_function=None
withnopot=None
m=None
def __init__(self,container,solver,elements,charge_operator,equation='potential',nonlocal_charge_potential_function=None):
"""
container: the container which describes the system.
solver: the solution of the system, e.g. generated by
my_container.lu_solver().
elements: the elements (from the system described by container) that
participate in the calculation. The property
fermi_energy_charge_dependence needs to be set.
charge_operator: The operator which gives the charge (without
\epsilon_0).
equation: there are three formulations of the problem: 'charge','potential' and 'charge_nonlocal'.
Default is 'potential'.
For 'charge', the element's charge_fermi_energy_dependence property is used.
For 'potential', the element's fermi_energy_charge_dependence property is used.
You don't have to set the property that you don't use.
For 'charge_nonlocal', you only have to provide the nonlocal_charge_potential_function
nonlocal_charge_potential_function: use together with equation='charge_nonlocal'.
Mind to calculate the volume charge density, not the area charge density!
The argument will be the difference vector Fermi_energy(x)-potential(x).
"""
self.container=container
self.solver=solver
self.elements=elements
self.charge_operator=charge_operator
if equation=='potential':
self.__equation=self.__functiontominimize
if equation=='charge':
self.__equation=self.__functiontominimize_charge
if equation=='charge_nonlocal':
self.__equation=self.__functiontominimize_charge_nonlocal
self.__nonlocal_charge_potential_function=nonlocal_charge_potential_function
def refresh_environment_contrib(self):
self.__reset_potential()
inhom=self.container.createinhomogeneity()
self.withnopot=self.container.charge(self.solver(inhom),self.charge_operator,self.elements)
def refresh_basisvecs(self,status=False):
self.__reset_potential()
if self.withnopot==None:
self.refresh_environment_contrib()
basisvecs=[]
for elem in self.elements:
if status:
print elem.index()
elem.potential=1
inhom=self.container.createinhomogeneity()
solution=self.solver(inhom)
basisvecs.append(self.container.charge(solution,self.charge_operator,self.elements)-self.withnopot)
elem.potential=0
self.m=numpy.array(basisvecs)
def __f(self,potvalues):
for x,elem in zip(potvalues,self.elements):
elem.potential=x
return self.__equation()
def __functiontominimize(self): #actually the root, not the minimum
potvec=numpy.array([elem.potential for elem in self.elements])
fermivec=numpy.array([elem.fermi_energy for elem in self.elements])
charges=numpy.dot(potvec,self.m)+self.withnopot
sqrtvec=[elem.fermi_energy_charge_dependence(charge)/Constants.elem_charge for elem,charge in zip(self.elements,charges)]
return -fermivec+potvec+sqrtvec
def __functiontominimize_charge(self):
"""
Minimize the charge difference instead of the potential difference.
"""
potvec=numpy.array([elem.potential for elem in self.elements])
fermivec=numpy.array([elem.fermi_energy for elem in self.elements])
charges_cl=numpy.dot(potvec,self.m)+self.withnopot
charges_qm=[elem.charge_fermi_energy_dependence(x*Constants.elem_charge) for x in fermivec-potvec]
"""
If rho is a functional:
charges_qm=charge_in_each_volume_element(fermivec-potvec)
Mind to calculate the volume charge density, not the area charge density!
"""
return charges_cl-charges_qm
def __functiontominimize_charge_nonlocal(self):
"""
Minimize the charge difference for a nonlocal connection between
potential and charge.
Works identical to __functiontominimize_charge, except for the
nonlocal qm charge calculation.
"""
potvec=numpy.array([elem.potential for elem in self.elements])
fermivec=numpy.array([elem.fermi_energy for elem in self.elements])
charges_cl=numpy.dot(potvec,self.m)+self.withnopot
charges_qm=self.__nonlocal_charge_potential_function(fermivec-potvec)
return charges_cl-charges_qm
def solve(self):
potvec=numpy.array([elem.potential for elem in self.elements])
loesung=scipy.optimize.root(self.__f,potvec,method='lm')
return loesung
def get_charge(self):
inhom=self.container.createinhomogeneity()
solution=self.solver(inhom)
charge=self.container.charge(solution,self.charge_operator,self.elements)
return charge
def get_potential(self):
potvec=numpy.array([elem.potential for elem in self.elements])
return potvec
def __reset_potential(self):
for elem in self.elements:
elem.potential=0
class FermiEnergyChargeDependence:
@staticmethod
[docs] def bulk_graphene(charge,grid_height=1e-9): #grid_height*volume density = area density
"""
How the Fermi energy of graphene depends on the charge density at 0K.
"""
return -Constants.v_fermi*Constants.hbar*math.copysign(math.sqrt(math.pi*math.fabs(charge*grid_height)/Constants.elem_charge),charge)
@staticmethod
def no_dependence(charge):
"""
How the Fermi energy of a "metal" depends on the charge density (it doesn't).
"""
return 0.
[docs]class BulkGrapheneWithTemperature:
"""
How the Fermi energy of graphene (Dirac cone) depends on the charge density at finite temperatures.
The algorithm includes a function inversion, which can either be done with interpolation (Ef_interp()) or
with a root search (Ef()). The interpolation is strongly suggested for higher speed.
If you use Q() often, you should probably also implement an interpolation
to improve speed.
The variables minEf, max Ef, Eftol, dEf are the parameters for the interpolation/root search. minEf and maxEf
have to span the relevant energy window around the equilibrium Fermi energy = 0. If steps occur, decrease Eftol/dEf.
Check Ef() if it's still working properly!!!
You will probably use Q() as charge_fermi_energy_dependence and
Ef() or Ef_interp() as fermi_energy_charge_dependence.
"""
T=None
grid_height=None
interp=None
minEf=-1.5e-18
maxEf=1.5e-18
Eftol=1e-25 #for root search
dEf=1e-22 #for interpolation
def __init__(self,T,grid_height,interpolation=True):
"""
If you don't want to use Ef_interp(), you can set
interpolation=False.
"""
self.T=T
self.grid_height=grid_height
Ef_grid=numpy.arange(self.maxEf,self.minEf,-self.dEf)
if interpolation:
Q_grid=[self.Q(i,T) for i in Ef_grid]
self.interp=scipy.interpolate.interp1d(Q_grid,Ef_grid)
def J1(self,eta):
return -float(sympy.mpmath.polylog(2,-sympy.exp(eta)))
def n(self,Ef,T):
return 2/math.pi*(Constants.k_B*T/(Constants.hbar*Constants.v_fermi))**2*self.J1(Ef/(Constants.k_B*T))
def p(self,Ef,T):
return 2/math.pi*(Constants.k_B*T/(Constants.hbar*Constants.v_fermi))**2*self.J1(-Ef/(Constants.k_B*T))
def Q(self,Ef,T=None):
"""
Note that Q is divided by grid_height to calculate a volume density
(instead of an area density).
"""
if T==None:
T=self.T
return Constants.elem_charge*(self.p(Ef,T)-self.n(Ef,T))/self.grid_height
def Ef(self,charge,T=None):
if T==None:
T=self.T
return scipy.optimize.brentq(lambda myEf: self.Q(myEf,T)-charge,self.minEf,self.maxEf,xtol=self.Eftol)
def Ef_interp(self,charge):
return self.interp(charge)