#Python source code. Requires [[:wikipedia:matplotlib|matplotlib]]. 
#Bonus features: 
1) Also makes entropy, energy, grand potential plots.
2) Also makes plots for Cl, and ionization of P, B impurities in silicon.

<syntaxhighlight lang="python">
"""
Plot various quantities related to thermal ionization of an atom,
calculated from simple model using grand canonical ensemble.
"""
from pylab import *
import matplotlib.transforms as transforms

plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.rc('font', serif=['Computer Modern'])
figtype = '.svg'
saveopts = {} #'bbox_inches':'tight'} #, 'transparent':True, 'frameon':True}
axsize = [0.2,0.16,0.79,0.77]

### Thermodynamic functions ###
def Omega(x, x_I, x_A, g_0, g_I, g_A):
    """
    Grand potential in terms of dimensionless parameters.
        x = -(\mu+e\phi)/kT, where
            \mu is chemical potential
            \phi is electrostatic potential of vacuum
            e is elementary charge
        x_I = E_I/kT, where E_I is ionization energy of atom
        x_A = E_A/kT, where E_A is electron affinity of atom
        g_0, g_I, g_A: degeneracies of neutral, oxidized, reduced states.
    
    Returns grand potential with neutral offsets left off, in units
    of kT:
        (Omega + mu N_0 - E_0)/kT
    """
    return -log(g_0 + g_I*exp(x-x_I) + g_A*exp(x_A-x))

def navg(x, x_I, x_A, g_0, g_I, g_A):
    """
    Average occupation number in terms of dimensionless parameters.
    (see Omega for parameters' meaning)
    
    Returns <N> - N_0
    
    The occupation number is given by <N> = - d\Omega/d\mu
    """
    return (-g_I*exp(x-x_I) + g_A*exp(x_A-x))/(g_0 + g_I*exp(x-x_I) + g_A*exp(x_A-x))

def entropy(x, x_I, x_A, g_0, g_I, g_A):
    """
    Entropy/k in terms of dimensionless parameters.    
    (see Omega for parameters' meaning)
    
    This function is calculated from the grand potential \Omega of this system,
    and is given by
        S/k = - d\Omega/d(kT)
    """
    t1 = log(g_0 + g_I*exp(x-x_I) + g_A*exp(x_A-x))
    t2 = -(g_I*(x-x_I)*exp(x-x_I) + g_A*(x_A-x)*exp(x_A-x))/(g_0 + g_I*exp(x-x_I) + g_A*exp(x_A-x))
    return t1+t2

### Figure maker template ###
def makefigs(name, W, kT = 1.,
    N_0 = 1, E_I = 1., E_A = 1.,
    g_0 = 1, g_I = 1, g_A = 1,
    bandgap=None):

    def makefig():
        fig = figure()
        fig.set_size_inches(3,3)
        fig.patch.set_alpha(0)
        ax = axes(axsize)
        xlim(amin(W),amax(W))
        trans = ax.get_xaxis_transform()
        if amin(W) <= E_A <= amax(W):
            axvline(E_A, color='gray', linestyle='dotted')
            text(E_A,1.01,r'$\Delta E_{\rm A}$',ha='center',va='bottom',transform=trans)
        if amin(W) <= E_I <= amax(W):
            axvline(E_I, color='gray', linestyle='dotted')
            text(E_I,1.01,r'$\Delta E_{\rm I}$',ha='center',va='bottom',transform=trans)
        if bandgap is None:
            # free atom terminology
            xlabel(r'$W = [-\mu-e\phi]$~(eV)')
        else:
            # semiconductor terminology
            xlabel(r'$\epsilon_{\rm C}-\mu$~(eV)')
            axvline(bandgap, color='gray', linestyle='dashed')
            text(bandgap,1.01,r'$\Delta E_{\rm gap}$',ha='center',va='bottom',transform=trans)
        return fig,ax

    N = navg(W/kT,E_I/kT,E_A/kT,g_0,g_I,g_A)
    S = entropy(W/kT,E_I/kT,E_A/kT,g_0,g_I,g_A)
    Om = Omega(W/kT,E_I/kT,E_A/kT,g_0,g_I,g_A)
    aveE = (Om + S)*kT # This is <E>-E_0, for the case when \mu = 0

    fig,ax = makefig()
    plot(W, N_0 + N, linewidth=1.5)
    Nmin, Nmax = N_0 + amin(N) - 0.2, N_0 + amax(N) + 0.2
    ylim(Nmin, Nmax)
    axhline(N_0, color='green', linestyle='solid', linewidth=0.5)
    ax.yaxis.set_ticks([t for t in range(N_0+10) if Nmin <= t <= Nmax])
    ylabel(r'$\langle N \rangle = -\frac{d\Omega}{d\mu}$')
    savefig('ionize_'+name+'_navg'+figtype,**saveopts)

    fig,ax = makefig()
    ylabel(r'$S/k = -\frac{d\Omega}{d(kT)}$')
    plot(W, S, linewidth=1.5)
    ylim(-0.1,1.7)
    savefig('ionize_'+name+'_entropy'+figtype,**saveopts)

    fig,ax = makefig()
    plot(W, Om*kT, linewidth=1.5)
    ax.autoscale(False)
    plot(W, -E_A + W - kT*log(g_A), color='gray', linewidth=0.5)
    plot(W, W*0 - kT*log(g_0), color='green', linewidth=0.5)
    plot(W, E_I - W - kT*log(g_I), color='gray', linewidth=0.5)
    text(0.5,0.95,r'(for $\mu = 0$)',ha='center',va='top',transform=ax.transAxes)
    ax.yaxis.set_label_coords(-0.17,0.5)
    ylabel(r'$\Omega - E_0$ (eV)')
    savefig('ionize_'+name+'_grand'+figtype,**saveopts)

    fig,ax = makefig()
    plot(W, aveE, color='b', linewidth=1.5)
    ax.autoscale(False)
    plot(W, -E_A + W, color='gray', linewidth=0.5)
    plot(W, W*0, color='green', linewidth=0.5)
    plot(W, E_I - W, color='gray', linewidth=0.5)
    text(0.5,0.95,r'(for $\mu = 0$)',ha='center',va='top',transform=ax.transAxes)
    ylabel(r'$\langle E \rangle - E_0$~(eV)')
    ax.yaxis.set_label_coords(-0.17,0.5)
    savefig('ionize_'+name+'_energy'+figtype,**saveopts)

### Specific data ###
makefigs('Cs', # free Cesium
    linspace(-0.4,5.2,541),
    kT = 8.61733238e-5 * 1500, #eV, 1500 K
    N_0 = 55,
    E_I = 3.89390, #eV, from WP:Ionization_energies_of_the_elements_(data_page)
    E_A = 0.47164, #eV, from WP:Electron_affinity_(data_page)
    g_0 = 2, # unpaired 6s electron spin degeneracy
    g_I = 1, # filled shells
    g_A = 1, # filled shells
    )

makefigs('Cl', # free Chlorine
    linspace(-0.4,5.2,541),
    kT = 8.61733238e-5 * 1500, #eV, 1500 K
    N_0 = 17,
    E_I = 12.96764, #eV, from WP:Ionization_energies_of_the_elements_(data_page)
    E_A = 3.612724, #eV, from WP:Electron_affinity_(data_page)
    g_0 = 2, # unpaired 3p hole spin degeneracy
    g_I = 1, # irrelevant placeholder value (no visible effect)
    g_A = 1, # filled shells
    )

# Below we try some ionization of dopants in silicon.
# The real behaviour is a bit more complicated than indicated here but this gives
# the conventional textbook model of dopant ionization.
#   See "Theory of shallow acceptor states in Si and Ge" by Schechter (1962)
#  also "The electronic structure of impurities and other point defects in semiconductors" by Pantelides (1978).

makefigs('Si-P', # Phosphorus in silicon (dopant)
    linspace(-0.1,1.25,261),
    kT = 8.61733238e-5 * 295, #eV, 295 K
    N_0 = 15,
    E_I = 0.045, #eV, from web
    E_A = -10, #eV, random large value to prevent ionization
    g_0 = 2, # 3sp^5 electron spin degeneracy, S=1/2 in this case
    g_I = 1, # half-filled shell of 3sp electrons... apparently nonmagnetic
    g_A = 1, # irrelevant placeholder value (no visible effect)
    bandgap = 1.1,
    )

makefigs('Si-B', # Boron in silicon (dopant)
    linspace(-0.1,1.25,261),
    kT = 8.61733238e-5 * 295, #eV, 295 K
    N_0 = 5,
    E_I = 10, #eV, random large value to prevent ionization
    E_A = 1.1-0.045, #eV, from web
    g_0 = 4, # 3sp^3 hole spin degeneracy: two possible orbital states (from two valence bands), each with S=1/2
    g_I = 1, # irrelevant placeholder value (no visible effect)
    g_A = 1, # half-filled shell of 3sp electrons... apparently nonmagnetic
    bandgap = 1.1,
    )