Source code for envtb.quantumcapacitance.quantumcapacitance

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