moscap.py

"""Stage 3: 1D MOS capacitor from scratch — the doorway to the MOSFET.

Structure: ideal metal gate / SiO2 (tox = 4.15 nm, sky130 1.8V-class) /
p-type silicon (NA = 1e17) / bulk contact. Equilibrium in silicon (no current
flows through the oxide), so only Poisson's equation is needed; the oxide
enters as a flux (Robin) boundary condition at the surface:

    eps_ox * (VG' - psi_s) / tox  =  eps_si * E_surface

Quasi-static (low-frequency) C-V by differentiating gate charge: C = dQg/dVG.

Checks against closed-form theory:
  - accumulation and strong-inversion capacitance -> Cox = eps_ox/tox
  - threshold voltage: VG at surface band bending 2*phiF matches
    VT = VFB + 2*phiF + sqrt(2*q*eps_si*NA*2*phiF)/Cox   (VFB = 0, ideal gate)
  - C-V minimum vs the high-frequency Cmin = series(Cox, eps_si/Wdm)
"""
from pathlib import Path

import numpy as np
from scipy.linalg import solve_banded
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

OUT = Path(__file__).parent

q = 1.602176634e-19
Vt = 0.025851997
ni = 1.0e10
eps_si = 11.7 * 8.85418782e-14
eps_ox = 3.9 * 8.85418782e-14

# --- the DESIGN ---------------------------------------------------------
TOX = 4.15e-7        # 4.15 nm gate oxide (sky130 1.8V device class)
NA = 1.0e17          # p-type body doping
DEPTH = 0.5e-4       # 0.5 um of silicon
VFB = 0.0            # ideal gate: no workfunction difference, no oxide charge

Cox = eps_ox / TOX
phiF = Vt * np.log(NA / ni)
psi_b = -Vt * np.arcsinh(NA / (2 * ni))          # bulk potential (intrinsic ref)

# geometric mesh: 0.1 nm at the surface (resolves the inversion layer),
# growing toward the bulk
h0, ratio = 1e-8, 1.035
hs = [h0]
while sum(hs) < DEPTH:
    hs.append(min(hs[-1] * ratio, 2e-6))
h = np.array(hs)
x = np.concatenate([[0.0], np.cumsum(h)])
N = x.size


def solve_gate(VG, psi):
    """Newton solve of Poisson in silicon with oxide Robin BC at node 0."""
    vg_eff = VG - VFB + psi_b          # gate electrode potential, intrinsic ref
    for _ in range(200):
        n = ni * np.exp(np.clip(psi / Vt, -60, 60))
        p = ni * np.exp(np.clip(-psi / Vt, -60, 60))
        rho = q * (p - n - NA)
        F = np.zeros(N)
        dF = np.zeros(N)
        # interior nodes (finite volume)
        F[1:-1] = eps_si * ((psi[2:] - psi[1:-1]) / h[1:] -
                            (psi[1:-1] - psi[:-2]) / h[:-1]) \
            + rho[1:-1] * (h[1:] + h[:-1]) / 2
        dF[1:-1] = -eps_si * (1 / h[1:] + 1 / h[:-1]) \
            - q * (n[1:-1] + p[1:-1]) / Vt * (h[1:] + h[:-1]) / 2
        # surface node: oxide flux + half control volume
        F[0] = Cox * (vg_eff - psi[0]) - eps_si * (psi[0] - psi[1]) / h[0] \
            + rho[0] * h[0] / 2
        dF[0] = -Cox - eps_si / h[0] - q * (n[0] + p[0]) / Vt * h[0] / 2
        # bulk node: Dirichlet
        F[-1] = 0.0
        dF[-1] = 1.0
        lower = np.zeros(N); upper = np.zeros(N)
        lower[1:-1] = eps_si / h[:-1]
        upper[1:-1] = eps_si / h[1:]
        upper[0] = eps_si / h[0]
        ab = np.zeros((3, N))
        ab[0, 1:] = upper[:-1]
        ab[1, :] = dF
        ab[2, :-1] = lower[1:]
        dpsi = solve_banded((1, 1), ab, -F)
        dpsi = np.clip(dpsi, -0.5, 0.5)
        psi = psi + dpsi
        if np.max(np.abs(dpsi)) < 1e-12:
            break
    return psi


# --- gate sweep ---------------------------------------------------------
VGs = np.arange(-3.0, 3.0001, 0.02)
psi = np.full(N, psi_b)
Qg, psis, snapshots = [], [], {}
for VG in VGs:
    psi = solve_gate(VG, psi)
    Qg.append(Cox * (VG - VFB + psi_b - psi[0]))     # charge on the gate
    psis.append(psi[0])
    for probe in (-2.0, 0.5, 2.0):
        if abs(VG - probe) < 1e-9:
            snapshots[probe] = psi.copy()

Qg, psis = np.array(Qg), np.array(psis)
C = np.gradient(Qg, VGs)                              # quasi-static C-V

# --- checks -------------------------------------------------------------
phis = psis - psi_b                                   # surface band bending
VT_analytic = VFB + 2 * phiF + np.sqrt(2 * q * eps_si * NA * 2 * phiF) / Cox
VT_sim = np.interp(2 * phiF, phis, VGs)               # VG where phis = 2 phiF
Wdm = np.sqrt(2 * eps_si * 2 * phiF / (q * NA))
Cmin_hf = Cox * (eps_si / Wdm) / (Cox + eps_si / Wdm)

C_acc = C[np.argmin(np.abs(VGs + 3.0))]
C_inv = C[np.argmin(np.abs(VGs - 3.0))]
print(f"Cox               : {Cox*1e9:.2f} nF/cm^2 (tox = {TOX*1e7:.2f} nm)")
print(f"C(acc, -3V)/Cox   : {C_acc/Cox:.3f}   (expect ->1 asymptotically)")
print(f"C(inv, +3V)/Cox   : {C_inv/Cox:.3f}   (expect ->1, low-frequency)")
print(f"min C / Cox       : {C.min()/Cox:.3f} vs analytic HF Cmin/Cox {Cmin_hf/Cox:.3f}")
print(f"threshold voltage : sim {VT_sim:.3f} V vs analytic {VT_analytic:.3f} V "
      f"(diff {1e3*abs(VT_sim-VT_analytic):.0f} mV)")
assert abs(C_acc / Cox - 1) < 0.05
assert abs(C_inv / Cox - 1) < 0.05
assert abs(VT_sim - VT_analytic) < 0.05
assert abs(C.min() - Cmin_hf) / Cmin_hf < 0.20

# --- plots --------------------------------------------------------------
fig, axes = plt.subplots(2, 2, figsize=(11, 7.5))

ax = axes[0, 0]
ax.plot(VGs, C / Cox, lw=2)
ax.axhline(Cmin_hf / Cox, color="gray", ls=":", label="analytic HF $C_{min}$")
ax.axvline(VT_analytic, color="tab:red", ls="--", lw=1,
           label=f"analytic $V_T$ = {VT_analytic:.2f} V")
ax.text(-1.9, 0.9, "accumulation"); ax.text(0.0, 0.45, "depletion")
ax.text(1.25, 0.9, "inversion")
ax.set_xlabel("$V_G$ (V)"); ax.set_ylabel("C / $C_{ox}$")
ax.set_title("quasi-static C-V"); ax.legend(loc="lower left"); ax.grid(alpha=0.3)

ax = axes[0, 1]
ax.plot(VGs, phis, lw=2)
ax.axhline(2 * phiF, color="tab:red", ls="--", lw=1, label="$2\\phi_F$")
ax.axvline(VT_sim, color="tab:red", ls=":", lw=1,
           label=f"sim $V_T$ = {VT_sim:.2f} V")
ax.set_xlabel("$V_G$ (V)"); ax.set_ylabel("surface band bending $\\phi_s$ (V)")
ax.set_title("band bending pins near $2\\phi_F$ past threshold")
ax.legend(); ax.grid(alpha=0.3)

ax = axes[1, 0]
for VG, color in ((-2.0, "tab:green"), (0.5, "tab:orange"), (2.0, "tab:blue")):
    psiV = snapshots[VG]
    nV = ni * np.exp(np.clip(psiV / Vt, -60, 60))
    pV = ni * np.exp(np.clip(-psiV / Vt, -60, 60))
    ax.semilogy(x * 1e7, nV, color=color, label=f"n, $V_G$={VG:g}")
    ax.semilogy(x * 1e7, pV, color=color, ls="--")
ax.axhline(NA, color="k", lw=0.6)
ax.set_xlim(0, 100)
ax.set_ylim(1e0, 1e20)
ax.set_xlabel("depth (nm)"); ax.set_ylabel("cm$^{-3}$")
ax.set_title("carriers: electrons (solid) / holes (dashed)")
ax.legend(fontsize=8); ax.grid(alpha=0.3)

ax = axes[1, 1]
ax.plot(VGs, -Qg / q * 1e-12, lw=2)
ax.set_xlabel("$V_G$ (V)")
ax.set_ylabel("semiconductor charge ($10^{12}$ q/cm$^2$)")
ax.set_title("total induced charge vs gate voltage")
ax.grid(alpha=0.3)

fig.suptitle("MOS capacitor from scratch — tox 4.15 nm, $N_A$=1e17 (sky130-class)",
             fontweight="bold")
fig.tight_layout()
fig.savefig(OUT / "moscap_results.png", dpi=130)
print(f"plot written: {OUT / 'moscap_results.png'}")
print("ALL CHECKS PASSED")