diff --git a/.gitignore b/.gitignore index 4655678..5e3c069 100644 --- a/.gitignore +++ b/.gitignore @@ -15,6 +15,7 @@ last_run_outputs/ *.last_log devsim-dev/ output_*/ +devices/*/output_*/ *.pkl *.txt nohup.out diff --git a/Makefile b/Makefile index 1626658..422ce04 100644 --- a/Makefile +++ b/Makefile @@ -10,9 +10,17 @@ export OPENBLAS_NUM_THREADS = 4 # Default simulation control options avalanche ?= false +btbt ?= false refine ?= false refine_v_step ?= 50.0 temp ?= 300.0 +pcad ?= false + +# Subdirectory device and output management +dev ?= Triac_rp +out ?= output_260625_01 +DEV_DIR = devices/$(dev) +OUT_DIR = $(DEV_DIR)/$(out) PYTHON := .venv/bin/python @@ -40,8 +48,7 @@ help: @echo " make show-conv - Print last few convergence step error details" @echo "" @echo "Output & Backup Rules:" - @echo " output_this_run/ - Current logs, checkpoints & visualization plots" - @echo " output_last_run/ - Previous run archive (triggered by new make mesh or sweep)" + @echo " devices/\$$(dev)/\$$(out)/ - Current logs, checkpoints & visualization plots" @echo " * Note: Rebuilding mesh via make mesh auto-archives output_this_run to prevent" @echo " loading old checkpoints on the updated grid structure." @echo "" @@ -54,8 +61,12 @@ help-detail: @echo "Detailed Command Parameters and Variables:" @echo "=============================================================================" @echo "Command Variables:" + @echo " dev= - Target device directory name (default: Triac_rp)" + @echo " out= - Target output folder name (default: output_260625_01)" @echo " avalanche=true|false - Toggle impact ionization (avalanche) model" @echo " Default: false (normal sweep without avalanche)" + @echo " btbt=true|false - Toggle band-to-band tunneling (BTBT) model" + @echo " Default: false (normal sweep without BTBT)" @echo " refine=true|false - Toggle dynamic adaptive refinement during sweep" @echo " Default: true (enable grid splitting at milestones)" @echo " refine_v_step= - Set voltage interval (V) to trigger dynamic refinement" @@ -66,11 +77,10 @@ help-detail: @echo " Default: automatically searches for latest checkpoints" @echo "" @echo "Usage Examples:" - @echo " make sweep avalanche=true" + @echo " make sweep dev=Triac_rp out=output_260625_01 pcad=true" + @echo " make sweep dev=Triac_rp avalanche=true" @echo " make sweep temp=350.0" - @echo " make sweep refine=false" - @echo " make sweep refine_v_step=25.0" - @echo " make resume checkpoint=output_this_run/seed_500V.pkl temp=350.0" + @echo " make resume checkpoint=devices/Triac_rp/output_260625_01/seed_500V.pkl temp=350.0" @echo " make resume-bg avalanche=true refine_v_step=30.0 temp=350.0" @echo "=============================================================================" @@ -79,55 +89,59 @@ help-detail: # 2. 執行基礎網格生成 # 3. 執行 run_refinement_2d.py 讀取基礎網格，求解電場並寫出新的 device_bgmesh.pos # 4. 再次執行 generate_mesh_2d.py，此時會自動載入 bgmesh 並輸出最終優化網格 device_2d.msh -mesh: device_config.py generate_mesh_2d.py generate_analytical_bgmesh.py - @echo ">>> [Mesh] 開始進行自適應網格重構流程..." - rm -f device_bgmesh.pos - $(PYTHON) generate_mesh_2d.py - $(PYTHON) generate_analytical_bgmesh.py - $(PYTHON) generate_mesh_2d.py - @echo ">>> [Mesh] 自適應優化網格生成完畢！(Saved: device_2d.msh)" +mesh: $(DEV_DIR)/device_config.py generate_mesh_2d.py generate_analytical_bgmesh.py + @echo ">>> [Mesh] 開始進行自適應網格重構流程 (dev=$(dev))..." + rm -f $(DEV_DIR)/device_bgmesh.pos + DEV_DIR=$(DEV_DIR) OUT_DIR=$(OUT_DIR)/ USE_PCAD=$(pcad) $(PYTHON) generate_mesh_2d.py + DEV_DIR=$(DEV_DIR) OUT_DIR=$(OUT_DIR)/ USE_PCAD=$(pcad) $(PYTHON) generate_analytical_bgmesh.py + DEV_DIR=$(DEV_DIR) OUT_DIR=$(OUT_DIR)/ USE_PCAD=$(pcad) $(PYTHON) generate_mesh_2d.py + @echo ">>> [Mesh] 自適應優化網格生成完畢！(Saved: $(DEV_DIR)/device_2d.msh)" # --- 熱平衡電位求解 --- # 依賴於對應的網格與求解腳本 -static: device_2d.msh solve_static_2d.py - @echo ">>> [Static] 求解零偏壓熱平衡狀態 (temp=$(temp))..." - TEMP=$(temp) $(PYTHON) solve_static_2d.py +static: $(DEV_DIR)/device_2d.msh solve_static_2d.py + @echo ">>> [Static] 求解零偏壓熱平衡狀態 (dev=$(dev), temp=$(temp))..." + DEV_DIR=$(DEV_DIR) OUT_DIR=$(OUT_DIR)/ USE_PCAD=$(pcad) TEMP=$(temp) $(PYTHON) solve_static_2d.py # --- 高壓偏壓掃描 --- # 依賴於對應的網格與掃描腳本 -sweep: device_2d.msh solve_sweep_recon.py - @echo ">>> [Sweep] 開始高壓偏壓漂移-擴散模擬 (avalanche=$(avalanche), refine=$(refine), refine_v_step=$(refine_v_step), temp=$(temp))..." - AVALANCHE=$(avalanche) REFINE=$(refine) REFINE_V_STEP=$(refine_v_step) TEMP=$(temp) $(PYTHON) solve_sweep_recon.py > sweeping.log 2>&1 +sweep: $(DEV_DIR)/device_2d.msh solve_sweep_recon.py + @echo ">>> [Sweep] 開始高壓偏壓漂移-擴散模擬 (dev=$(dev), out=$(out), avalanche=$(avalanche), btbt=$(btbt), refine=$(refine), refine_v_step=$(refine_v_step), temp=$(temp))..." + mkdir -p $(OUT_DIR) + DEV_DIR=$(DEV_DIR) OUT_DIR=$(OUT_DIR)/ USE_PCAD=$(pcad) AVALANCHE=$(avalanche) BTBT=$(btbt) REFINE=$(refine) REFINE_V_STEP=$(refine_v_step) TEMP=$(temp) $(PYTHON) -u solve_sweep_recon.py > $(OUT_DIR)/sweeping.log 2>&1 resume: - @echo ">>> [Resume] 從指定的 Checkpoint ($(checkpoint)) 或最新的自動備份接續掃描 (avalanche=$(avalanche), refine=$(refine), refine_v_step=$(refine_v_step), temp=$(temp))..." - AVALANCHE=$(avalanche) REFINE=$(refine) REFINE_V_STEP=$(refine_v_step) TEMP=$(temp) $(PYTHON) resume_run.py $(checkpoint) >> sweeping.log 2>&1 + @echo ">>> [Resume] 從指定的 Checkpoint ($(checkpoint)) 或最新的自動備份接續掃描 (dev=$(dev), out=$(out), avalanche=$(avalanche), btbt=$(btbt), refine=$(refine), refine_v_step=$(refine_v_step), temp=$(temp))..." + mkdir -p $(OUT_DIR) + DEV_DIR=$(DEV_DIR) OUT_DIR=$(OUT_DIR)/ USE_PCAD=$(pcad) AVALANCHE=$(avalanche) BTBT=$(btbt) REFINE=$(refine) REFINE_V_STEP=$(refine_v_step) TEMP=$(temp) $(PYTHON) -u resume_run.py $(checkpoint) >> $(OUT_DIR)/sweeping.log 2>&1 resume-bg: - @echo ">>> [Resume-BG] 在背景從指定的 Checkpoint ($(checkpoint)) 或最新的自動備份接續掃描 (avalanche=$(avalanche), refine=$(refine), refine_v_step=$(refine_v_step), temp=$(temp))..." - AVALANCHE=$(avalanche) REFINE=$(refine) REFINE_V_STEP=$(refine_v_step) TEMP=$(temp) nohup $(PYTHON) resume_run.py $(checkpoint) >> sweeping.log 2>&1 & + @echo ">>> [Resume-BG] 在背景從指定的 Checkpoint ($(checkpoint)) 或最新的自動備份接續掃描 (dev=$(dev), out=$(out), avalanche=$(avalanche), btbt=$(btbt), refine=$(refine), refine_v_step=$(refine_v_step), temp=$(temp))..." + mkdir -p $(OUT_DIR) + DEV_DIR=$(DEV_DIR) OUT_DIR=$(OUT_DIR)/ USE_PCAD=$(pcad) AVALANCHE=$(avalanche) BTBT=$(btbt) REFINE=$(refine) REFINE_V_STEP=$(refine_v_step) TEMP=$(temp) nohup $(PYTHON) -u resume_run.py $(checkpoint) >> $(OUT_DIR)/sweeping.log 2>&1 & # --- 萃取與監控收斂曲線 --- show-conv: - @if [ -f sweeping.log ]; then \ - awk '/Iteration:/ {printf "Iteration %s：", $$2} /Device:/ {print $$4}' sweeping.log | tail -n 10; \ + @if [ -f $(OUT_DIR)/sweeping.log ]; then \ + awk '/Iteration:/ {printf "Iteration %s：", $$2} /Device:/ {print $$4}' $(OUT_DIR)/sweeping.log | tail -n 10; \ else \ - echo "sweeping.log does not exist."; \ + echo "$(OUT_DIR)/sweeping.log does not exist."; \ fi monitor: - @if [ -f sweeping.log ]; then \ - tail -f sweeping.log | awk '/Iteration:/ {printf "Iteration %s：", $$2; fflush()} /Device:/ {print $$4; fflush()}'; \ + @if [ -f $(OUT_DIR)/sweeping.log ]; then \ + tail -f $(OUT_DIR)/sweeping.log | awk '/Iteration:/ {printf "Iteration %s：", $$2; fflush()} /Device:/ {print $$4; fflush()}'; \ else \ - echo "sweeping.log does not exist."; \ + echo "$(OUT_DIR)/sweeping.log does not exist."; \ fi # --- 網格依賴規則 --- # 當沒有 device_2d.msh 或 device_config.py 有更動時，自動觸發 mesh 流程 -device_2d.msh: device_config.py generate_mesh_2d.py generate_analytical_bgmesh.py - $(MAKE) mesh +$(DEV_DIR)/device_2d.msh: $(DEV_DIR)/device_config.py generate_mesh_2d.py generate_analytical_bgmesh.py + $(MAKE) mesh dev=$(dev) clean: @echo ">>> 清除暫存與網格檔案..." rm -f *.msh *.pos *.tec *.png *.csv *.vtm *.vtu *.visit rm -rf __pycache__ physics/__pycache__ + diff --git a/device_pcad_config.py b/device_pcad_config.py new file mode 100644 index 0000000..02f295c --- /dev/null +++ b/device_pcad_config.py @@ -0,0 +1,340 @@ +# device_pcad_config.py +# Process-Step Oriented 2D Doping Configuration File +# This file coexists with device_config.py and keeps the original doping path intact. + +import os +import sys +import numpy as np +import math + +# All units in cm (1 um = 1e-4 cm) +um = 1e-4 + +# Import geometric and simulation parameters from the original device_config +from device_config import ( + W_DEVICE, H_SI, T_OX, H_MOLD, W_SIDE_MOLD, W_SIM, + VIA_WIDTH, VIA_P11_X, VIA_P13_X, + MT1_FP1_X1, MT1_FP1_X2, MT1_FP2_X1, MT1_FP2_X2, + SIM_NAME +) + +# Define mask openings (Horizontal coordinate ranges, positive X-axis. Mirroring is handled dynamically) +masks = { + 'p_well_mask': [(75.0 * um, 100.0 * um)], + 'p12_mask': [(120.0 * um, 130.0 * um)], + 'p13_mask': [(150.0 * um, 255.0 * um)], + 'nplus_mask': [(164.0 * um, 185.0 * um)], + 'mring_mask': [(340.0 * um, W_DEVICE)], +} + +# --- Process Database --- +# Implant Range & Straggle in Silicon (Energy in keV, Rp and dRp in microns) +IMPLANT_DB = { + 'Boron': { + 'energy': np.array([10.0, 20.0, 30.0, 50.0, 80.0, 100.0, 150.0, 200.0]), + 'Rp': np.array([0.04, 0.07, 0.10, 0.16, 0.24, 0.30, 0.42, 0.55]), + 'dRp': np.array([0.015, 0.025, 0.035, 0.050, 0.063, 0.070, 0.085, 0.100]) + }, + 'Phosphorus': { + 'energy': np.array([10.0, 20.0, 30.0, 50.0, 80.0, 100.0, 150.0, 200.0]), + 'Rp': np.array([0.015, 0.028, 0.040, 0.065, 0.100, 0.120, 0.180, 0.240]), + 'dRp': np.array([0.007, 0.012, 0.016, 0.025, 0.038, 0.045, 0.060, 0.075]) + }, + 'Arsenic': { + 'energy': np.array([10.0, 20.0, 30.0, 50.0, 80.0, 100.0, 150.0, 200.0]), + 'Rp': np.array([0.010, 0.018, 0.025, 0.038, 0.058, 0.070, 0.100, 0.130]), + 'dRp': np.array([0.004, 0.007, 0.009, 0.014, 0.020, 0.025, 0.035, 0.045]) + } +} + +# Diffusion Arrhenius parameters in Silicon (D0 in cm^2/s, Ea in eV) +DIFFUSION_DB = { + 'Boron': {'D0': 1.0, 'Ea': 3.46}, + 'Phosphorus': {'D0': 10.5, 'Ea': 3.69}, + 'Arsenic': {'D0': 0.32, 'Ea': 3.56} +} + +def get_implant_params(dopant, energy_kev): + """Interpolate Rp and dRp from database for a given energy in keV. Returns values in cm.""" + db = IMPLANT_DB.get(dopant, IMPLANT_DB['Boron']) + e = np.clip(energy_kev, db['energy'][0], db['energy'][-1]) + Rp_um = np.interp(e, db['energy'], db['Rp']) + dRp_um = np.interp(e, db['energy'], db['dRp']) + return Rp_um * 1e-4, dRp_um * 1e-4 + +def get_diffusion_coefficient(dopant, temp_c): + """Calculate the diffusion coefficient D in cm^2/s at temperature temp_c in Celsius.""" + db = DIFFUSION_DB.get(dopant, DIFFUSION_DB['Boron']) + T_k = temp_c + 273.15 + k_B = 8.617333262145e-5 # eV/K + D0 = db['D0'] + Ea = db['Ea'] + return D0 * np.exp(-Ea / (k_B * T_k)) + +# Define process steps for the 2D doping profile +process_steps = [ + # Step 0: Substrate Wafer + { + 'enabled': True, + 'name': 'Substrate', + 'method': 'Substrate', + 'type': 'n', + 'doping': 1.0e14, # N_SUB + }, + # Step 1: P11 Well Implant & Diffusion + { + 'enabled': True, + 'name': 'P11_Well', + 'method': 'Implant_and_Diffuse', + 'type': 'p', + 'dopant': 'Boron', + 'energy': 80.0, # keV + 'dose': 1.77245385091e12, # cm^-2 + 'mask': 'p_well_mask', + 'temp': 1150.0, # Celsius + 'time': 465.54227, # minutes + 'lateral_ratio': 0.8, + }, + # Step 2: P12 Well Implant & Diffusion + { + 'enabled': True, + 'name': 'P12_Well', + 'method': 'Implant_and_Diffuse', + 'type': 'p', + 'dopant': 'Boron', + 'energy': 80.0, + 'dose': 6.6467019409e11, + 'mask': 'p12_mask', + 'temp': 1150.0, + 'time': 465.54227, + 'lateral_ratio': 0.8, + }, + # Step 3: P13 Well Implant & Diffusion + { + 'enabled': True, + 'name': 'P13_Well', + 'method': 'Implant_and_Diffuse', + 'type': 'p', + 'dopant': 'Boron', + 'energy': 80.0, + 'dose': 1.77245385091e12, + 'mask': 'p13_mask', + 'temp': 1150.0, + 'time': 465.54227, + 'lateral_ratio': 0.8, + }, + # Step 4: N+ Source Active Region + { + 'enabled': True, + 'name': 'NPlus_Active', + 'method': 'Implant_and_Diffuse', + 'type': 'n', + 'dopant': 'Phosphorus', + 'energy': 50.0, + 'dose': 2.65868077636e11, + 'mask': 'nplus_mask', + 'temp': 1000.0, + 'time': 34.108586, + 'lateral_ratio': 0.666666666667, + }, + # Step 5: MRING Guard Ring Active Region + { + 'enabled': True, + 'name': 'MRing_Active', + 'method': 'Implant_and_Diffuse', + 'type': 'n', + 'dopant': 'Phosphorus', + 'energy': 50.0, + 'dose': 2.65868077636e11, + 'mask': 'mring_mask', + 'temp': 1000.0, + 'time': 34.108586, + 'lateral_ratio': 0.666666666667, + } +] + +def resolve_steps(): + """ + Resolves the physical process steps into analytical model parameters + (peak, vdiff, hdiff, x_ranges) for profile evaluation. + """ + resolved = [] + + # 1. Resolve substrate first + sub_step = [s for s in process_steps if s['method'] == 'Substrate'][0] + if sub_step.get('enabled', True): + resolved.append({ + 'name': sub_step.get('name', 'Substrate'), + 'method': 'Substrate', + 'type': sub_step['type'], + 'doping': sub_step['doping'] + }) + + # 2. Resolve implant/diffusion steps + for step in process_steps: + if step['method'] == 'Substrate' or not step.get('enabled', True): + continue + + name = step.get('name') + type_ = step['type'] + dopant = step['dopant'] + method = step['method'] + + # Look up mask coordinate ranges + mask_name = step['mask'] + x_ranges = masks.get(mask_name, []) + + # Get implant properties + energy = step.get('energy', 80.0) + dose = step.get('dose', 1.0e12) + Rp, dRp = get_implant_params(dopant, energy) + + # Get thermal diffusion properties + temp = step.get('temp', 1000.0) + time_min = step.get('time', 30.0) + + D = get_diffusion_coefficient(dopant, temp) + time_sec = time_min * 60.0 + Dt = D * time_sec + + # Standard deviation incorporating both implant straggle and thermal budget diffusion + vdiff = math.sqrt(2.0 * (dRp ** 2) + 4.0 * Dt) + + # Horizontal standard deviation based on lateral ratio + lateral_ratio = step.get('lateral_ratio', 0.8) + hdiff = vdiff * lateral_ratio + + # Calculate peak concentration matching the integrated dose for erfc profile: + # peak = dose / (vdiff * (sqrt(pi) / 2)) + peak = dose / (vdiff * (math.sqrt(math.pi) / 2.0)) + + resolved.append({ + 'name': name, + 'method': 'Implant_and_Diffuse', + 'type': type_, + 'peak': peak, + 'vdiff': vdiff, + 'hdiff': hdiff, + 'x_ranges': x_ranges + }) + + return resolved + +def apply_pcad_doping_2d(device, region="Silicon"): + """ + Parses the process steps list and registers the corresponding + 2D analytical doping models in DEVSIM. + """ + import devsim + + donor_terms = [] + acceptor_terms = [] + + steps = resolve_steps() + + # 1. Setup Substrate baseline first + sub_step = [s for s in steps if s['method'] == 'Substrate'][0] + sub_doping = sub_step['doping'] + sub_type = sub_step['type'] + + devsim.node_model(device=device, region=region, name="nD_sub", equation=f"{sub_doping}") + + if sub_type == 'n': + donor_terms.append("nD_sub") + else: + acceptor_terms.append("nD_sub") + + # 2. Iterate and build expressions for Implant/Diffusion steps + idx = 1 + for step in steps: + if step['method'] == 'Substrate': + continue + + name = step['name'] + type_ = step['type'] + peak = step['peak'] + x_ranges = step['x_ranges'] + vdiff = step['vdiff'] + hdiff = step['hdiff'] + + # Construct analytical 2D erfc expression for each window range, including mirroring + expr_terms = [] + for x1, x2 in x_ranges: + # Right side (positive X) + expr_terms.append(f"erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))") + # Left side (negative X, mirrored: x1 -> -x1 and x2 -> -x2, so window is [-x2, -x1]) + expr_terms.append(f"erfc(y / {vdiff}) * 0.5 * (erf((x - ({-x2})) / {hdiff}) - erf((x - ({-x1})) / {hdiff}))") + + # Combine terms and multiply by peak concentration + combined_expr = f"{peak} * (" + " + ".join(expr_terms) + ")" + + model_name = f"nA_{name}" if type_ == 'p' else f"nD_{name}" + devsim.node_model(device=device, region=region, name=model_name, equation=combined_expr) + + if type_ == 'n': + donor_terms.append(model_name) + else: + acceptor_terms.append(model_name) + + idx += 1 + + # 3. Combine separate models into global Donors and Acceptors fields in DEVSIM + donor_equation = " + ".join(donor_terms) + acceptor_equation = "1e10 + " + " + ".join(acceptor_terms) + + devsim.node_model(device=device, region=region, name="Donors", equation=donor_equation) + devsim.node_model(device=device, region=region, name="Acceptors", equation=acceptor_equation) + devsim.node_model(device=device, region=region, name="NetDoping", equation="Donors - Acceptors") + devsim.node_model(device=device, region=region, name="LogNetDoping", equation="asinh(NetDoping / 2.0) / log(10.0)") + +def get_pcad_doping_val_2d(x, y): + """ + Evaluates the 2D process-step doping profile analytically on numpy arrays x and y. + """ + import numpy as np + import math + + erf_vec = np.vectorize(math.erf) + erfc_vec = np.vectorize(math.erfc) + + donors = np.zeros_like(x, dtype=float) + acceptors = np.zeros_like(x, dtype=float) + + steps = resolve_steps() + + # 1. Setup Substrate baseline first + sub_step = [s for s in steps if s['method'] == 'Substrate'][0] + sub_doping = sub_step['doping'] + sub_type = sub_step['type'] + if sub_type == 'n': + donors += sub_doping + else: + acceptors += sub_doping + + # 2. Iterate and build expressions for Implant/Diffusion steps + for step in steps: + if step['method'] == 'Substrate': + continue + + type_ = step['type'] + peak = step['peak'] + x_ranges = step['x_ranges'] + vdiff = step['vdiff'] + hdiff = step['hdiff'] + + prof = np.zeros_like(x, dtype=float) + for x1, x2 in x_ranges: + # Right side + prof += erfc_vec(y / vdiff) * 0.5 * (erf_vec((x - x1) / hdiff) - erf_vec((x - x2) / hdiff)) + # Left side (mirrored) + prof += erfc_vec(y / vdiff) * 0.5 * (erf_vec((x - (-x2)) / hdiff) - erf_vec((x - (-x1)) / hdiff)) + + if type_ == 'n': + donors += peak * prof + else: + acceptors += peak * prof + + # Acceptor base offset 1e10 + acceptors += 1e10 + + return donors - acceptors diff --git a/devices/LDMOS/device_config.py b/devices/LDMOS/device_config.py new file mode 100644 index 0000000..5b753f7 --- /dev/null +++ b/devices/LDMOS/device_config.py @@ -0,0 +1,69 @@ +# device_config.py +# All units in cm (1 um = 1e-4 cm) +um = 1e-4 + +# --- Geometric Dimensions --- +W_DEVICE = 400.0 * um # Half-width of the device (400 x 2 total width) +H_SI = 200.0 * um # Silicon substrate thickness +T_OX = 2.0 * um # Oxide thickness +H_MOLD = 20.0 * um # Molding compound thickness (above oxide) +W_SIDE_MOLD = 20.0 * um # Molding compound width on the sides +W_SIM = W_DEVICE + W_SIDE_MOLD # Half-width of the total simulation domain + +# --- P-well parameters (p11, p12, p13) --- +P_WELL_DEPTH = 5.0 * um # 5 um depth for all P-wells + +# P-well X boundaries (Right half, will be mirrored for left half) +P11_X1 = 75.0 * um +P11_X2 = 100.0 * um + +P12_X1 = 120.0 * um +P12_X2 = 130.0 * um + +P13_X1 = 150.0 * um +P13_X2 = 255.0 * um + +# --- N+ region parameters --- +NPLUS_DEPTH = 0.1 * um # 1 um depth for all N+ regions + +# N+ X boundaries (Right half, mirrored for left half) +NPLUS_X1 = 164.0 * um +NPLUS_X2 = 185.0 * um + +# MRING X boundaries (Right half, mirrored for left half) +MRING_X1 = 340.0 * um +MRING_X2 = 356.0 * um + +# --- Doping Concentrations (cm^-3) --- +N_SUB = 1.0e14 # 70 ~ 90 ohm cm (5.5e13) +P11_PEAK = 8.0e15 +P12_PEAK = 3.0e15 +P13_PEAK = 8.0e15 +NPLUS_PEAK = 2.0e16 + +# --- Doping Gradient / Diffusion Widths --- +# P-well gradient widths +P_WELL_VDDIFF = 2.5 * um # Vertical gradient width (characteristic depth) +P_WELL_HDDIFF = 2.0 * um # Horizontal (lateral) gradient width + +# N+ gradient widths +NPLUS_VDDIFF = 0.15 * um # Vertical gradient width +NPLUS_HDDIFF = 0.1 * um # Horizontal (lateral) gradient width + +# --- Contact Vias Width and Positions (Right half, mirrored for left) --- +VIA_WIDTH = 3.0 * um + +# Contact via center positions +VIA_P11_X = 87.5 * um +VIA_P13_X = 174.5 * um + +# --- Metal Field Plate X boundaries (Right half, mirrored for left) --- +MT1_FP1_X1 = 30.0 * um +MT1_FP1_X2 = 186.0 * um + +MT1_FP2_X1 = 250.0 * um +MT1_FP2_X2 = 295.0 * um + +# --- Simulation Metadata --- +SIM_NAME = "symmetric LDMOS 20260625" + diff --git a/devices/LDMOS/device_pcad_config.py b/devices/LDMOS/device_pcad_config.py new file mode 100644 index 0000000..02f295c --- /dev/null +++ b/devices/LDMOS/device_pcad_config.py @@ -0,0 +1,340 @@ +# device_pcad_config.py +# Process-Step Oriented 2D Doping Configuration File +# This file coexists with device_config.py and keeps the original doping path intact. + +import os +import sys +import numpy as np +import math + +# All units in cm (1 um = 1e-4 cm) +um = 1e-4 + +# Import geometric and simulation parameters from the original device_config +from device_config import ( + W_DEVICE, H_SI, T_OX, H_MOLD, W_SIDE_MOLD, W_SIM, + VIA_WIDTH, VIA_P11_X, VIA_P13_X, + MT1_FP1_X1, MT1_FP1_X2, MT1_FP2_X1, MT1_FP2_X2, + SIM_NAME +) + +# Define mask openings (Horizontal coordinate ranges, positive X-axis. Mirroring is handled dynamically) +masks = { + 'p_well_mask': [(75.0 * um, 100.0 * um)], + 'p12_mask': [(120.0 * um, 130.0 * um)], + 'p13_mask': [(150.0 * um, 255.0 * um)], + 'nplus_mask': [(164.0 * um, 185.0 * um)], + 'mring_mask': [(340.0 * um, W_DEVICE)], +} + +# --- Process Database --- +# Implant Range & Straggle in Silicon (Energy in keV, Rp and dRp in microns) +IMPLANT_DB = { + 'Boron': { + 'energy': np.array([10.0, 20.0, 30.0, 50.0, 80.0, 100.0, 150.0, 200.0]), + 'Rp': np.array([0.04, 0.07, 0.10, 0.16, 0.24, 0.30, 0.42, 0.55]), + 'dRp': np.array([0.015, 0.025, 0.035, 0.050, 0.063, 0.070, 0.085, 0.100]) + }, + 'Phosphorus': { + 'energy': np.array([10.0, 20.0, 30.0, 50.0, 80.0, 100.0, 150.0, 200.0]), + 'Rp': np.array([0.015, 0.028, 0.040, 0.065, 0.100, 0.120, 0.180, 0.240]), + 'dRp': np.array([0.007, 0.012, 0.016, 0.025, 0.038, 0.045, 0.060, 0.075]) + }, + 'Arsenic': { + 'energy': np.array([10.0, 20.0, 30.0, 50.0, 80.0, 100.0, 150.0, 200.0]), + 'Rp': np.array([0.010, 0.018, 0.025, 0.038, 0.058, 0.070, 0.100, 0.130]), + 'dRp': np.array([0.004, 0.007, 0.009, 0.014, 0.020, 0.025, 0.035, 0.045]) + } +} + +# Diffusion Arrhenius parameters in Silicon (D0 in cm^2/s, Ea in eV) +DIFFUSION_DB = { + 'Boron': {'D0': 1.0, 'Ea': 3.46}, + 'Phosphorus': {'D0': 10.5, 'Ea': 3.69}, + 'Arsenic': {'D0': 0.32, 'Ea': 3.56} +} + +def get_implant_params(dopant, energy_kev): + """Interpolate Rp and dRp from database for a given energy in keV. Returns values in cm.""" + db = IMPLANT_DB.get(dopant, IMPLANT_DB['Boron']) + e = np.clip(energy_kev, db['energy'][0], db['energy'][-1]) + Rp_um = np.interp(e, db['energy'], db['Rp']) + dRp_um = np.interp(e, db['energy'], db['dRp']) + return Rp_um * 1e-4, dRp_um * 1e-4 + +def get_diffusion_coefficient(dopant, temp_c): + """Calculate the diffusion coefficient D in cm^2/s at temperature temp_c in Celsius.""" + db = DIFFUSION_DB.get(dopant, DIFFUSION_DB['Boron']) + T_k = temp_c + 273.15 + k_B = 8.617333262145e-5 # eV/K + D0 = db['D0'] + Ea = db['Ea'] + return D0 * np.exp(-Ea / (k_B * T_k)) + +# Define process steps for the 2D doping profile +process_steps = [ + # Step 0: Substrate Wafer + { + 'enabled': True, + 'name': 'Substrate', + 'method': 'Substrate', + 'type': 'n', + 'doping': 1.0e14, # N_SUB + }, + # Step 1: P11 Well Implant & Diffusion + { + 'enabled': True, + 'name': 'P11_Well', + 'method': 'Implant_and_Diffuse', + 'type': 'p', + 'dopant': 'Boron', + 'energy': 80.0, # keV + 'dose': 1.77245385091e12, # cm^-2 + 'mask': 'p_well_mask', + 'temp': 1150.0, # Celsius + 'time': 465.54227, # minutes + 'lateral_ratio': 0.8, + }, + # Step 2: P12 Well Implant & Diffusion + { + 'enabled': True, + 'name': 'P12_Well', + 'method': 'Implant_and_Diffuse', + 'type': 'p', + 'dopant': 'Boron', + 'energy': 80.0, + 'dose': 6.6467019409e11, + 'mask': 'p12_mask', + 'temp': 1150.0, + 'time': 465.54227, + 'lateral_ratio': 0.8, + }, + # Step 3: P13 Well Implant & Diffusion + { + 'enabled': True, + 'name': 'P13_Well', + 'method': 'Implant_and_Diffuse', + 'type': 'p', + 'dopant': 'Boron', + 'energy': 80.0, + 'dose': 1.77245385091e12, + 'mask': 'p13_mask', + 'temp': 1150.0, + 'time': 465.54227, + 'lateral_ratio': 0.8, + }, + # Step 4: N+ Source Active Region + { + 'enabled': True, + 'name': 'NPlus_Active', + 'method': 'Implant_and_Diffuse', + 'type': 'n', + 'dopant': 'Phosphorus', + 'energy': 50.0, + 'dose': 2.65868077636e11, + 'mask': 'nplus_mask', + 'temp': 1000.0, + 'time': 34.108586, + 'lateral_ratio': 0.666666666667, + }, + # Step 5: MRING Guard Ring Active Region + { + 'enabled': True, + 'name': 'MRing_Active', + 'method': 'Implant_and_Diffuse', + 'type': 'n', + 'dopant': 'Phosphorus', + 'energy': 50.0, + 'dose': 2.65868077636e11, + 'mask': 'mring_mask', + 'temp': 1000.0, + 'time': 34.108586, + 'lateral_ratio': 0.666666666667, + } +] + +def resolve_steps(): + """ + Resolves the physical process steps into analytical model parameters + (peak, vdiff, hdiff, x_ranges) for profile evaluation. + """ + resolved = [] + + # 1. Resolve substrate first + sub_step = [s for s in process_steps if s['method'] == 'Substrate'][0] + if sub_step.get('enabled', True): + resolved.append({ + 'name': sub_step.get('name', 'Substrate'), + 'method': 'Substrate', + 'type': sub_step['type'], + 'doping': sub_step['doping'] + }) + + # 2. Resolve implant/diffusion steps + for step in process_steps: + if step['method'] == 'Substrate' or not step.get('enabled', True): + continue + + name = step.get('name') + type_ = step['type'] + dopant = step['dopant'] + method = step['method'] + + # Look up mask coordinate ranges + mask_name = step['mask'] + x_ranges = masks.get(mask_name, []) + + # Get implant properties + energy = step.get('energy', 80.0) + dose = step.get('dose', 1.0e12) + Rp, dRp = get_implant_params(dopant, energy) + + # Get thermal diffusion properties + temp = step.get('temp', 1000.0) + time_min = step.get('time', 30.0) + + D = get_diffusion_coefficient(dopant, temp) + time_sec = time_min * 60.0 + Dt = D * time_sec + + # Standard deviation incorporating both implant straggle and thermal budget diffusion + vdiff = math.sqrt(2.0 * (dRp ** 2) + 4.0 * Dt) + + # Horizontal standard deviation based on lateral ratio + lateral_ratio = step.get('lateral_ratio', 0.8) + hdiff = vdiff * lateral_ratio + + # Calculate peak concentration matching the integrated dose for erfc profile: + # peak = dose / (vdiff * (sqrt(pi) / 2)) + peak = dose / (vdiff * (math.sqrt(math.pi) / 2.0)) + + resolved.append({ + 'name': name, + 'method': 'Implant_and_Diffuse', + 'type': type_, + 'peak': peak, + 'vdiff': vdiff, + 'hdiff': hdiff, + 'x_ranges': x_ranges + }) + + return resolved + +def apply_pcad_doping_2d(device, region="Silicon"): + """ + Parses the process steps list and registers the corresponding + 2D analytical doping models in DEVSIM. + """ + import devsim + + donor_terms = [] + acceptor_terms = [] + + steps = resolve_steps() + + # 1. Setup Substrate baseline first + sub_step = [s for s in steps if s['method'] == 'Substrate'][0] + sub_doping = sub_step['doping'] + sub_type = sub_step['type'] + + devsim.node_model(device=device, region=region, name="nD_sub", equation=f"{sub_doping}") + + if sub_type == 'n': + donor_terms.append("nD_sub") + else: + acceptor_terms.append("nD_sub") + + # 2. Iterate and build expressions for Implant/Diffusion steps + idx = 1 + for step in steps: + if step['method'] == 'Substrate': + continue + + name = step['name'] + type_ = step['type'] + peak = step['peak'] + x_ranges = step['x_ranges'] + vdiff = step['vdiff'] + hdiff = step['hdiff'] + + # Construct analytical 2D erfc expression for each window range, including mirroring + expr_terms = [] + for x1, x2 in x_ranges: + # Right side (positive X) + expr_terms.append(f"erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))") + # Left side (negative X, mirrored: x1 -> -x1 and x2 -> -x2, so window is [-x2, -x1]) + expr_terms.append(f"erfc(y / {vdiff}) * 0.5 * (erf((x - ({-x2})) / {hdiff}) - erf((x - ({-x1})) / {hdiff}))") + + # Combine terms and multiply by peak concentration + combined_expr = f"{peak} * (" + " + ".join(expr_terms) + ")" + + model_name = f"nA_{name}" if type_ == 'p' else f"nD_{name}" + devsim.node_model(device=device, region=region, name=model_name, equation=combined_expr) + + if type_ == 'n': + donor_terms.append(model_name) + else: + acceptor_terms.append(model_name) + + idx += 1 + + # 3. Combine separate models into global Donors and Acceptors fields in DEVSIM + donor_equation = " + ".join(donor_terms) + acceptor_equation = "1e10 + " + " + ".join(acceptor_terms) + + devsim.node_model(device=device, region=region, name="Donors", equation=donor_equation) + devsim.node_model(device=device, region=region, name="Acceptors", equation=acceptor_equation) + devsim.node_model(device=device, region=region, name="NetDoping", equation="Donors - Acceptors") + devsim.node_model(device=device, region=region, name="LogNetDoping", equation="asinh(NetDoping / 2.0) / log(10.0)") + +def get_pcad_doping_val_2d(x, y): + """ + Evaluates the 2D process-step doping profile analytically on numpy arrays x and y. + """ + import numpy as np + import math + + erf_vec = np.vectorize(math.erf) + erfc_vec = np.vectorize(math.erfc) + + donors = np.zeros_like(x, dtype=float) + acceptors = np.zeros_like(x, dtype=float) + + steps = resolve_steps() + + # 1. Setup Substrate baseline first + sub_step = [s for s in steps if s['method'] == 'Substrate'][0] + sub_doping = sub_step['doping'] + sub_type = sub_step['type'] + if sub_type == 'n': + donors += sub_doping + else: + acceptors += sub_doping + + # 2. Iterate and build expressions for Implant/Diffusion steps + for step in steps: + if step['method'] == 'Substrate': + continue + + type_ = step['type'] + peak = step['peak'] + x_ranges = step['x_ranges'] + vdiff = step['vdiff'] + hdiff = step['hdiff'] + + prof = np.zeros_like(x, dtype=float) + for x1, x2 in x_ranges: + # Right side + prof += erfc_vec(y / vdiff) * 0.5 * (erf_vec((x - x1) / hdiff) - erf_vec((x - x2) / hdiff)) + # Left side (mirrored) + prof += erfc_vec(y / vdiff) * 0.5 * (erf_vec((x - (-x2)) / hdiff) - erf_vec((x - (-x1)) / hdiff)) + + if type_ == 'n': + donors += peak * prof + else: + acceptors += peak * prof + + # Acceptor base offset 1e10 + acceptors += 1e10 + + return donors - acceptors diff --git a/devices/Triac_rp/device_config.py b/devices/Triac_rp/device_config.py new file mode 100644 index 0000000..74f493e --- /dev/null +++ b/devices/Triac_rp/device_config.py @@ -0,0 +1,69 @@ +# device_config.py +# All units in cm (1 um = 1e-4 cm) +um = 1e-4 + +# --- Geometric Dimensions --- +W_DEVICE = 356.0 * um # Half-width of the device (356 x 2 total width) +H_SI = 200.0 * um # Silicon substrate thickness +T_OX = 2.0 * um # Oxide thickness +H_MOLD = 100.0 * um # Molding compound thickness (above oxide) +W_SIDE_MOLD = 100.0 * um # Molding compound width on the sides +W_SIM = W_DEVICE + W_SIDE_MOLD # Half-width of the total simulation domain + +# --- P-well parameters (p11, p12, p13) --- +P_WELL_DEPTH = 5.0 * um # 5 um depth for all P-wells + +# P-well X boundaries (Right half, will be mirrored for left half) +P11_X1 = 75.0 * um +P11_X2 = 100.0 * um + +P12_X1 = 120.0 * um +P12_X2 = 130.0 * um + +P13_X1 = 150.0 * um +P13_X2 = 255.0 * um + +# --- N+ region parameters --- +NPLUS_DEPTH = 0.1 * um # 1 um depth for all N+ regions + +# N+ X boundaries (Right half, mirrored for left half) +NPLUS_X1 = 164.0 * um +NPLUS_X2 = 185.0 * um + +# MRING X boundaries (Right half, mirrored for left half) +MRING_X1 = 340.0 * um +MRING_X2 = 356.0 * um + +# --- Doping Concentrations (cm^-3) --- +N_SUB = 1.0e14 # 70 ~ 90 ohm cm (5.5e13) +P11_PEAK = 8.0e15 +P12_PEAK = 3.0e15 +P13_PEAK = 8.0e15 +NPLUS_PEAK = 2.0e16 + +# --- Doping Gradient / Diffusion Widths --- +# P-well gradient widths +P_WELL_VDDIFF = 2.5 * um # Vertical gradient width (characteristic depth) +P_WELL_HDDIFF = 2.0 * um # Horizontal (lateral) gradient width + +# N+ gradient widths +NPLUS_VDDIFF = 0.15 * um # Vertical gradient width +NPLUS_HDDIFF = 0.1 * um # Horizontal (lateral) gradient width + +# --- Contact Vias Width and Positions (Right half, mirrored for left) --- +VIA_WIDTH = 10.0 * um + +# Contact via center positions +VIA_P11_X = 87.5 * um +VIA_P13_X = 174.5 * um + +# --- Metal Field Plate X boundaries (Right half, mirrored for left) --- +MT1_FP1_X1 = 30.0 * um +MT1_FP1_X2 = 186.0 * um + +MT1_FP2_X1 = 250.0 * um +MT1_FP2_X2 = 295.0 * um + +# --- Simulation Metadata --- +SIM_NAME = "p-doping 8&3e15 20260616" + diff --git a/devices/Triac_rp/device_pcad_config.py b/devices/Triac_rp/device_pcad_config.py new file mode 100644 index 0000000..02f295c --- /dev/null +++ b/devices/Triac_rp/device_pcad_config.py @@ -0,0 +1,340 @@ +# device_pcad_config.py +# Process-Step Oriented 2D Doping Configuration File +# This file coexists with device_config.py and keeps the original doping path intact. + +import os +import sys +import numpy as np +import math + +# All units in cm (1 um = 1e-4 cm) +um = 1e-4 + +# Import geometric and simulation parameters from the original device_config +from device_config import ( + W_DEVICE, H_SI, T_OX, H_MOLD, W_SIDE_MOLD, W_SIM, + VIA_WIDTH, VIA_P11_X, VIA_P13_X, + MT1_FP1_X1, MT1_FP1_X2, MT1_FP2_X1, MT1_FP2_X2, + SIM_NAME +) + +# Define mask openings (Horizontal coordinate ranges, positive X-axis. Mirroring is handled dynamically) +masks = { + 'p_well_mask': [(75.0 * um, 100.0 * um)], + 'p12_mask': [(120.0 * um, 130.0 * um)], + 'p13_mask': [(150.0 * um, 255.0 * um)], + 'nplus_mask': [(164.0 * um, 185.0 * um)], + 'mring_mask': [(340.0 * um, W_DEVICE)], +} + +# --- Process Database --- +# Implant Range & Straggle in Silicon (Energy in keV, Rp and dRp in microns) +IMPLANT_DB = { + 'Boron': { + 'energy': np.array([10.0, 20.0, 30.0, 50.0, 80.0, 100.0, 150.0, 200.0]), + 'Rp': np.array([0.04, 0.07, 0.10, 0.16, 0.24, 0.30, 0.42, 0.55]), + 'dRp': np.array([0.015, 0.025, 0.035, 0.050, 0.063, 0.070, 0.085, 0.100]) + }, + 'Phosphorus': { + 'energy': np.array([10.0, 20.0, 30.0, 50.0, 80.0, 100.0, 150.0, 200.0]), + 'Rp': np.array([0.015, 0.028, 0.040, 0.065, 0.100, 0.120, 0.180, 0.240]), + 'dRp': np.array([0.007, 0.012, 0.016, 0.025, 0.038, 0.045, 0.060, 0.075]) + }, + 'Arsenic': { + 'energy': np.array([10.0, 20.0, 30.0, 50.0, 80.0, 100.0, 150.0, 200.0]), + 'Rp': np.array([0.010, 0.018, 0.025, 0.038, 0.058, 0.070, 0.100, 0.130]), + 'dRp': np.array([0.004, 0.007, 0.009, 0.014, 0.020, 0.025, 0.035, 0.045]) + } +} + +# Diffusion Arrhenius parameters in Silicon (D0 in cm^2/s, Ea in eV) +DIFFUSION_DB = { + 'Boron': {'D0': 1.0, 'Ea': 3.46}, + 'Phosphorus': {'D0': 10.5, 'Ea': 3.69}, + 'Arsenic': {'D0': 0.32, 'Ea': 3.56} +} + +def get_implant_params(dopant, energy_kev): + """Interpolate Rp and dRp from database for a given energy in keV. Returns values in cm.""" + db = IMPLANT_DB.get(dopant, IMPLANT_DB['Boron']) + e = np.clip(energy_kev, db['energy'][0], db['energy'][-1]) + Rp_um = np.interp(e, db['energy'], db['Rp']) + dRp_um = np.interp(e, db['energy'], db['dRp']) + return Rp_um * 1e-4, dRp_um * 1e-4 + +def get_diffusion_coefficient(dopant, temp_c): + """Calculate the diffusion coefficient D in cm^2/s at temperature temp_c in Celsius.""" + db = DIFFUSION_DB.get(dopant, DIFFUSION_DB['Boron']) + T_k = temp_c + 273.15 + k_B = 8.617333262145e-5 # eV/K + D0 = db['D0'] + Ea = db['Ea'] + return D0 * np.exp(-Ea / (k_B * T_k)) + +# Define process steps for the 2D doping profile +process_steps = [ + # Step 0: Substrate Wafer + { + 'enabled': True, + 'name': 'Substrate', + 'method': 'Substrate', + 'type': 'n', + 'doping': 1.0e14, # N_SUB + }, + # Step 1: P11 Well Implant & Diffusion + { + 'enabled': True, + 'name': 'P11_Well', + 'method': 'Implant_and_Diffuse', + 'type': 'p', + 'dopant': 'Boron', + 'energy': 80.0, # keV + 'dose': 1.77245385091e12, # cm^-2 + 'mask': 'p_well_mask', + 'temp': 1150.0, # Celsius + 'time': 465.54227, # minutes + 'lateral_ratio': 0.8, + }, + # Step 2: P12 Well Implant & Diffusion + { + 'enabled': True, + 'name': 'P12_Well', + 'method': 'Implant_and_Diffuse', + 'type': 'p', + 'dopant': 'Boron', + 'energy': 80.0, + 'dose': 6.6467019409e11, + 'mask': 'p12_mask', + 'temp': 1150.0, + 'time': 465.54227, + 'lateral_ratio': 0.8, + }, + # Step 3: P13 Well Implant & Diffusion + { + 'enabled': True, + 'name': 'P13_Well', + 'method': 'Implant_and_Diffuse', + 'type': 'p', + 'dopant': 'Boron', + 'energy': 80.0, + 'dose': 1.77245385091e12, + 'mask': 'p13_mask', + 'temp': 1150.0, + 'time': 465.54227, + 'lateral_ratio': 0.8, + }, + # Step 4: N+ Source Active Region + { + 'enabled': True, + 'name': 'NPlus_Active', + 'method': 'Implant_and_Diffuse', + 'type': 'n', + 'dopant': 'Phosphorus', + 'energy': 50.0, + 'dose': 2.65868077636e11, + 'mask': 'nplus_mask', + 'temp': 1000.0, + 'time': 34.108586, + 'lateral_ratio': 0.666666666667, + }, + # Step 5: MRING Guard Ring Active Region + { + 'enabled': True, + 'name': 'MRing_Active', + 'method': 'Implant_and_Diffuse', + 'type': 'n', + 'dopant': 'Phosphorus', + 'energy': 50.0, + 'dose': 2.65868077636e11, + 'mask': 'mring_mask', + 'temp': 1000.0, + 'time': 34.108586, + 'lateral_ratio': 0.666666666667, + } +] + +def resolve_steps(): + """ + Resolves the physical process steps into analytical model parameters + (peak, vdiff, hdiff, x_ranges) for profile evaluation. + """ + resolved = [] + + # 1. Resolve substrate first + sub_step = [s for s in process_steps if s['method'] == 'Substrate'][0] + if sub_step.get('enabled', True): + resolved.append({ + 'name': sub_step.get('name', 'Substrate'), + 'method': 'Substrate', + 'type': sub_step['type'], + 'doping': sub_step['doping'] + }) + + # 2. Resolve implant/diffusion steps + for step in process_steps: + if step['method'] == 'Substrate' or not step.get('enabled', True): + continue + + name = step.get('name') + type_ = step['type'] + dopant = step['dopant'] + method = step['method'] + + # Look up mask coordinate ranges + mask_name = step['mask'] + x_ranges = masks.get(mask_name, []) + + # Get implant properties + energy = step.get('energy', 80.0) + dose = step.get('dose', 1.0e12) + Rp, dRp = get_implant_params(dopant, energy) + + # Get thermal diffusion properties + temp = step.get('temp', 1000.0) + time_min = step.get('time', 30.0) + + D = get_diffusion_coefficient(dopant, temp) + time_sec = time_min * 60.0 + Dt = D * time_sec + + # Standard deviation incorporating both implant straggle and thermal budget diffusion + vdiff = math.sqrt(2.0 * (dRp ** 2) + 4.0 * Dt) + + # Horizontal standard deviation based on lateral ratio + lateral_ratio = step.get('lateral_ratio', 0.8) + hdiff = vdiff * lateral_ratio + + # Calculate peak concentration matching the integrated dose for erfc profile: + # peak = dose / (vdiff * (sqrt(pi) / 2)) + peak = dose / (vdiff * (math.sqrt(math.pi) / 2.0)) + + resolved.append({ + 'name': name, + 'method': 'Implant_and_Diffuse', + 'type': type_, + 'peak': peak, + 'vdiff': vdiff, + 'hdiff': hdiff, + 'x_ranges': x_ranges + }) + + return resolved + +def apply_pcad_doping_2d(device, region="Silicon"): + """ + Parses the process steps list and registers the corresponding + 2D analytical doping models in DEVSIM. + """ + import devsim + + donor_terms = [] + acceptor_terms = [] + + steps = resolve_steps() + + # 1. Setup Substrate baseline first + sub_step = [s for s in steps if s['method'] == 'Substrate'][0] + sub_doping = sub_step['doping'] + sub_type = sub_step['type'] + + devsim.node_model(device=device, region=region, name="nD_sub", equation=f"{sub_doping}") + + if sub_type == 'n': + donor_terms.append("nD_sub") + else: + acceptor_terms.append("nD_sub") + + # 2. Iterate and build expressions for Implant/Diffusion steps + idx = 1 + for step in steps: + if step['method'] == 'Substrate': + continue + + name = step['name'] + type_ = step['type'] + peak = step['peak'] + x_ranges = step['x_ranges'] + vdiff = step['vdiff'] + hdiff = step['hdiff'] + + # Construct analytical 2D erfc expression for each window range, including mirroring + expr_terms = [] + for x1, x2 in x_ranges: + # Right side (positive X) + expr_terms.append(f"erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))") + # Left side (negative X, mirrored: x1 -> -x1 and x2 -> -x2, so window is [-x2, -x1]) + expr_terms.append(f"erfc(y / {vdiff}) * 0.5 * (erf((x - ({-x2})) / {hdiff}) - erf((x - ({-x1})) / {hdiff}))") + + # Combine terms and multiply by peak concentration + combined_expr = f"{peak} * (" + " + ".join(expr_terms) + ")" + + model_name = f"nA_{name}" if type_ == 'p' else f"nD_{name}" + devsim.node_model(device=device, region=region, name=model_name, equation=combined_expr) + + if type_ == 'n': + donor_terms.append(model_name) + else: + acceptor_terms.append(model_name) + + idx += 1 + + # 3. Combine separate models into global Donors and Acceptors fields in DEVSIM + donor_equation = " + ".join(donor_terms) + acceptor_equation = "1e10 + " + " + ".join(acceptor_terms) + + devsim.node_model(device=device, region=region, name="Donors", equation=donor_equation) + devsim.node_model(device=device, region=region, name="Acceptors", equation=acceptor_equation) + devsim.node_model(device=device, region=region, name="NetDoping", equation="Donors - Acceptors") + devsim.node_model(device=device, region=region, name="LogNetDoping", equation="asinh(NetDoping / 2.0) / log(10.0)") + +def get_pcad_doping_val_2d(x, y): + """ + Evaluates the 2D process-step doping profile analytically on numpy arrays x and y. + """ + import numpy as np + import math + + erf_vec = np.vectorize(math.erf) + erfc_vec = np.vectorize(math.erfc) + + donors = np.zeros_like(x, dtype=float) + acceptors = np.zeros_like(x, dtype=float) + + steps = resolve_steps() + + # 1. Setup Substrate baseline first + sub_step = [s for s in steps if s['method'] == 'Substrate'][0] + sub_doping = sub_step['doping'] + sub_type = sub_step['type'] + if sub_type == 'n': + donors += sub_doping + else: + acceptors += sub_doping + + # 2. Iterate and build expressions for Implant/Diffusion steps + for step in steps: + if step['method'] == 'Substrate': + continue + + type_ = step['type'] + peak = step['peak'] + x_ranges = step['x_ranges'] + vdiff = step['vdiff'] + hdiff = step['hdiff'] + + prof = np.zeros_like(x, dtype=float) + for x1, x2 in x_ranges: + # Right side + prof += erfc_vec(y / vdiff) * 0.5 * (erf_vec((x - x1) / hdiff) - erf_vec((x - x2) / hdiff)) + # Left side (mirrored) + prof += erfc_vec(y / vdiff) * 0.5 * (erf_vec((x - (-x2)) / hdiff) - erf_vec((x - (-x1)) / hdiff)) + + if type_ == 'n': + donors += peak * prof + else: + acceptors += peak * prof + + # Acceptor base offset 1e10 + acceptors += 1e10 + + return donors - acceptors diff --git a/dynamic_refine.py b/dynamic_refine.py index a70dd58..046619c 100644 --- a/dynamic_refine.py +++ b/dynamic_refine.py @@ -11,7 +11,7 @@ from physics.new_physics import * import generate_mesh_2d OUT_DIR = "output_this_run/" -def setup_physics_for_device(device, is_avalanche_enabled=False): +def setup_physics_for_device(device, is_avalanche_enabled=False, is_btbt_enabled=False): def CreateOxidePotentialOnlyContact(device, region, contact): contact_bias = GetContactBiasName(contact) contact_model = f"Potential - {contact_bias}" @@ -205,16 +205,40 @@ def setup_physics_for_device(device, is_avalanche_enabled=False): # Avalanche generation model (enabled/disabled by refine_and_interpolate caller) if is_avalanche_enabled: CreateAvalancheGeneration(device, "Silicon", opts['Jn'], opts['Jp']) + # BTBT generation model + if is_btbt_enabled: + CreateBTBTGeneration(device, "Silicon") av_model_n = "AvalancheGeneration" if is_avalanche_enabled else "" av_model_p = "AvalancheGeneration_p" if is_avalanche_enabled else "" + btbt_model_n = "BTBTGeneration" if is_btbt_enabled else "" + btbt_model_p = "BTBTGeneration_p" if is_btbt_enabled else "" + + if av_model_n and btbt_model_n: + from physics.model_create import CreateEdgeModel, CreateEdgeModelDerivatives + CreateEdgeModel(device, "Silicon", "CombinedGeneration", "AvalancheGeneration + BTBTGeneration") + CreateEdgeModel(device, "Silicon", "CombinedGeneration_p", "AvalancheGeneration_p + BTBTGeneration_p") + for i in ("Potential", "Electrons", "Holes"): + CreateEdgeModelDerivatives(device, "Silicon", "CombinedGeneration", "AvalancheGeneration + BTBTGeneration", i) + CreateEdgeModelDerivatives(device, "Silicon", "CombinedGeneration_p", "AvalancheGeneration_p + BTBTGeneration_p", i) + gen_model_n = "CombinedGeneration" + gen_model_p = "CombinedGeneration_p" + elif av_model_n: + gen_model_n = av_model_n + gen_model_p = av_model_p + elif btbt_model_n: + gen_model_n = btbt_model_n + gen_model_p = btbt_model_p + else: + gen_model_n = "" + gen_model_p = "" # 預設以 positive (full Newton) 方式註冊連續方程式 # refine_and_interpolate 在插值後會臨時切換至 log_damp 做 Stage 1 預處理 devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons", - time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=av_model_n, + time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=gen_model_n, variable_update="positive", node_model="ElectronGeneration", min_error=1e5) devsim.equation(device=device, region="Silicon", name="HoleContinuityEquation", variable_name="Holes", - time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=av_model_p, + time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=gen_model_p, variable_update="positive", node_model="HoleGeneration", min_error=1e5) devsim.node_model(device=device, region="Silicon", name="LogElectrons", equation="log(Electrons + 1e-10) / log(10.0)") @@ -388,7 +412,7 @@ def enforce_contact_boundary_conditions(device_name): except Exception as ex: print(f" Error enforcing {reg} contact boundary conditions: {ex}") -def refine_and_interpolate(device_old, v_bias, is_avalanche_enabled=False, time_log=None, out_dir=None): +def refine_and_interpolate(device_old, v_bias, is_avalanche_enabled=False, is_btbt_enabled=False, time_log=None, out_dir=None): if out_dir is None: out_dir = OUT_DIR import time @@ -465,7 +489,7 @@ def refine_and_interpolate(device_old, v_bias, is_avalanche_enabled=False, time_ # 5. 載入新網格並設定物理與 solutions devsim.create_gmsh_mesh(mesh=device_new_name, file=mesh_out_path) - opts = setup_physics_for_device(device_new_name, is_avalanche_enabled=is_avalanche_enabled) + opts = setup_physics_for_device(device_new_name, is_avalanche_enabled=is_avalanche_enabled, is_btbt_enabled=is_btbt_enabled) # 6. Apply bias to contacts of the new device for c in ["MT1_Si", "MT1_P12_Si", "MT1_Ox", "MT1_Mold"]: @@ -764,14 +788,35 @@ def refine_and_interpolate(device_old, v_bias, is_avalanche_enabled=False, time_ # ========================================== # Stage 1: Fully-coupled log_damp # ========================================== - # Re-register Electron and Hole Continuity equations in Silicon and contacts av_model_n = "AvalancheGeneration" if is_avalanche_enabled else "" av_model_p = "AvalancheGeneration_p" if is_avalanche_enabled else "" + btbt_model_n = "BTBTGeneration" if is_btbt_enabled else "" + btbt_model_p = "BTBTGeneration_p" if is_btbt_enabled else "" + + if av_model_n and btbt_model_n: + from physics.model_create import CreateEdgeModel, CreateEdgeModelDerivatives + CreateEdgeModel(device_new_name, "Silicon", "CombinedGeneration", "AvalancheGeneration + BTBTGeneration") + CreateEdgeModel(device_new_name, "Silicon", "CombinedGeneration_p", "AvalancheGeneration_p + BTBTGeneration_p") + for i in ("Potential", "Electrons", "Holes"): + CreateEdgeModelDerivatives(device_new_name, "Silicon", "CombinedGeneration", "AvalancheGeneration + BTBTGeneration", i) + CreateEdgeModelDerivatives(device_new_name, "Silicon", "CombinedGeneration_p", "AvalancheGeneration_p + BTBTGeneration_p", i) + gen_model_n = "CombinedGeneration" + gen_model_p = "CombinedGeneration_p" + elif av_model_n: + gen_model_n = av_model_n + gen_model_p = av_model_p + elif btbt_model_n: + gen_model_n = btbt_model_n + gen_model_p = btbt_model_p + else: + gen_model_n = "" + gen_model_p = "" + devsim.equation(device=device_new_name, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons", - time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=av_model_n, + time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=gen_model_n, variable_update="log_damp", node_model="ElectronGeneration", min_error=1e5) devsim.equation(device=device_new_name, region="Silicon", name="HoleContinuityEquation", variable_name="Holes", - time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=av_model_p, + time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=gen_model_p, variable_update="log_damp", node_model="HoleGeneration", min_error=1e5) for c in ["MT1_Si", "MT2_Si", "MT1_P12_Si", "MT2_P12_Si"]: contact_electrons_name = f"{c}nodeelectrons" @@ -840,11 +885,33 @@ def refine_and_interpolate(device_old, v_bias, is_avalanche_enabled=False, time_ # ========================================== av_model_n = "AvalancheGeneration" if is_avalanche_enabled else "" av_model_p = "AvalancheGeneration_p" if is_avalanche_enabled else "" + btbt_model_n = "BTBTGeneration" if is_btbt_enabled else "" + btbt_model_p = "BTBTGeneration_p" if is_btbt_enabled else "" + + if av_model_n and btbt_model_n: + from physics.model_create import CreateEdgeModel, CreateEdgeModelDerivatives + CreateEdgeModel(device_new_name, "Silicon", "CombinedGeneration", "AvalancheGeneration + BTBTGeneration") + CreateEdgeModel(device_new_name, "Silicon", "CombinedGeneration_p", "AvalancheGeneration_p + BTBTGeneration_p") + for i in ("Potential", "Electrons", "Holes"): + CreateEdgeModelDerivatives(device_new_name, "Silicon", "CombinedGeneration", "AvalancheGeneration + BTBTGeneration", i) + CreateEdgeModelDerivatives(device_new_name, "Silicon", "CombinedGeneration_p", "AvalancheGeneration_p + BTBTGeneration_p", i) + gen_model_n = "CombinedGeneration" + gen_model_p = "CombinedGeneration_p" + elif av_model_n: + gen_model_n = av_model_n + gen_model_p = av_model_p + elif btbt_model_n: + gen_model_n = btbt_model_n + gen_model_p = btbt_model_p + else: + gen_model_n = "" + gen_model_p = "" + devsim.equation(device=device_new_name, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons", - time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=av_model_n, + time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=gen_model_n, variable_update="positive", node_model="ElectronGeneration", min_error=1e5) devsim.equation(device=device_new_name, region="Silicon", name="HoleContinuityEquation", variable_name="Holes", - time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=av_model_p, + time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=gen_model_p, variable_update="positive", node_model="HoleGeneration", min_error=1e5) # Restore PotentialEquation variable update to default with min_error=1e-3 devsim.equation(device=device_new_name, region="Silicon", name="PotentialEquation", variable_name="Potential", @@ -905,7 +972,7 @@ def refine_and_interpolate(device_old, v_bias, is_avalanche_enabled=False, time_ devsim.create_gmsh_mesh(mesh=device_old, file="device_2d.msh") except Exception: pass - setup_physics_for_device(device_old, is_avalanche_enabled=is_avalanche_enabled) + setup_physics_for_device(device_old, is_avalanche_enabled=is_avalanche_enabled, is_btbt_enabled=is_btbt_enabled) raise devsim.error(f"Precision Newton solve on refined mesh did not converge: {e}") print("Convergence on refined mesh achieved successfully!") diff --git a/generate_analytical_bgmesh.py b/generate_analytical_bgmesh.py index 25bfa4a..47bc7db 100644 --- a/generate_analytical_bgmesh.py +++ b/generate_analytical_bgmesh.py @@ -4,41 +4,48 @@ import math import sys import os -sys.path.append("/home/pchan/devsim2026") +DEV_DIR = os.environ.get("DEV_DIR", "devices/Triac_rp") +sys.path.insert(0, os.path.abspath(DEV_DIR)) from device_config import * # Vectorize math functions for fast numpy operations erf_vec = np.vectorize(math.erf) erfc_vec = np.vectorize(math.erfc) -def erfc_doping(x, y, peak, x1, x2, hdiff, vdiff): - return peak * erfc_vec(y / vdiff) * 0.5 * (erf_vec((x - x1) / hdiff) - erf_vec((x - x2) / hdiff)) +if os.environ.get("USE_PCAD", "false").lower() == "true": + from device_pcad_config import get_pcad_doping_val_2d + def get_doping_val(x, y): + return get_pcad_doping_val_2d(x, y) +else: + def erfc_doping(x, y, peak, x1, x2, hdiff, vdiff): + return peak * erfc_vec(y / vdiff) * 0.5 * (erf_vec((x - x1) / hdiff) - erf_vec((x - x2) / hdiff)) -def get_doping_val(x, y): - # Donors - nD = N_SUB - nD += erfc_doping(x, y, NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) - nD += erfc_doping(x, y, NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF) - nD += erfc_doping(x, y, NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) - nD += erfc_doping(x, y, NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF) - - # Acceptors - nA = 1e10 - nA += erfc_doping(x, y, P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) - nA += erfc_doping(x, y, P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) - nA += erfc_doping(x, y, P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) - nA += erfc_doping(x, y, P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) - nA += erfc_doping(x, y, P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) - nA += erfc_doping(x, y, P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) - - return nD - nA + def get_doping_val(x, y): + # Donors + nD = N_SUB + nD += erfc_doping(x, y, NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) + nD += erfc_doping(x, y, NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF) + nD += erfc_doping(x, y, NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) + nD += erfc_doping(x, y, NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF) + + # Acceptors + nA = 1e10 + nA += erfc_doping(x, y, P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + nA += erfc_doping(x, y, P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + nA += erfc_doping(x, y, P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + nA += erfc_doping(x, y, P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + nA += erfc_doping(x, y, P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + nA += erfc_doping(x, y, P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + + return nD - nA def generate_analytical_bgmesh(): device = "device_2d" # Load the mesh generated by first pass of generate_mesh_2d.py (coarse base mesh) - print("Loading base mesh: device_2d.msh...") - devsim.create_gmsh_mesh(mesh=device, file="device_2d.msh") + mesh_path = os.path.join(DEV_DIR, "device_2d.msh") + print(f"Loading base mesh: {mesh_path}...") + devsim.create_gmsh_mesh(mesh=device, file=mesh_path) devsim.add_gmsh_region(mesh=device, gmsh_name="Silicon", region="Silicon", material="Silicon") devsim.add_gmsh_region(mesh=device, gmsh_name="Oxide", region="Oxide", material="Oxide") devsim.add_gmsh_region(mesh=device, gmsh_name="Molding", region="Molding", material="Molding") @@ -51,7 +58,7 @@ def generate_analytical_bgmesh(): N_offset = 1.0e10 # Intrinsic concentration offset to avoid division by zero G_ref = 0.2 / um # Reference relative gradient (0.2 um^-1) for transition smoothness - bgmesh_path = "device_bgmesh.pos" + bgmesh_path = os.path.join(DEV_DIR, "device_bgmesh.pos") with open(bgmesh_path, "w") as f: f.write('View "background mesh" {\n') diff --git a/generate_mesh_2d.py b/generate_mesh_2d.py index 14a8a34..0f7fc2a 100644 --- a/generate_mesh_2d.py +++ b/generate_mesh_2d.py @@ -1,6 +1,9 @@ import gmsh import numpy as np import os +import sys +DEV_DIR = os.environ.get("DEV_DIR", "devices/Triac_rp") +sys.path.insert(0, os.path.abspath(DEV_DIR)) from device_config import * def create_mesh(y_box_max=12.0*um, y_medium_max=20.0*um, mesh_out="device_2d.msh", bgmesh_pos="device_bgmesh.pos"): @@ -348,4 +351,7 @@ def create_mesh(y_box_max=12.0*um, y_medium_max=20.0*um, mesh_out="device_2d.msh print("Mesh generation complete! Saved as device_2d.msh.") if __name__ == "__main__": - create_mesh() + DEV_DIR = os.environ.get("DEV_DIR", "devices/Triac_rp") + mesh_out = os.path.join(DEV_DIR, "device_2d.msh") + bgmesh_pos = os.path.join(DEV_DIR, "device_bgmesh.pos") + create_mesh(mesh_out=mesh_out, bgmesh_pos=bgmesh_pos) diff --git a/gui1d/app.py b/gui1d/app.py index e3fd839..5c75934 100644 --- a/gui1d/app.py +++ b/gui1d/app.py @@ -17,7 +17,7 @@ if ROOT_DIR not in sys.path: import subprocess import pickle -def build_and_solve_1d(bias_target, substrate_type, substrate_doping, length, process_steps, enable_avalanche, area_cm2): +def build_and_solve_1d(bias_target, substrate_type, substrate_doping, length, process_steps, enable_avalanche, enable_btbt, area_cm2): # Run the simulation in a separate python process to ensure thread-safety and prevent DEVSIM C++ segment faults temp_dir = os.path.join(ROOT_DIR, "gui1d", "temp_runs") os.makedirs(temp_dir, exist_ok=True) @@ -33,6 +33,7 @@ def build_and_solve_1d(bias_target, substrate_type, substrate_doping, length, pr 'length': length, 'process_steps': process_steps, 'enable_avalanche': enable_avalanche, + 'enable_btbt': enable_btbt, 'area_cm2': area_cm2 } @@ -49,11 +50,13 @@ def build_and_solve_1d(bias_target, substrate_type, substrate_doping, length, pr result = pickle.load(f) return result - except Exception as e: + except subprocess.CalledProcessError as e: err_msg = str(e) - if 'res_proc' in locals() and res_proc.stderr: - err_msg += "\n" + res_proc.stderr + if e.stderr: + err_msg += "\n" + e.stderr raise RuntimeError(f"Simulation process failed: {err_msg}") + except Exception as e: + raise RuntimeError(f"Simulation process failed: {str(e)}") finally: for fpath in (in_file, out_file): if os.path.exists(fpath): @@ -359,16 +362,21 @@ with st.sidebar: st.session_state[f"step_{idx}_cs"] = float(step.get('cs', 1e19)) st.session_state[f"step_{idx}_temp"] = float(step.get('temp', 1000.0)) st.session_state[f"step_{idx}_time"] = float(step.get('time', 60.0)) + st.session_state[f"step_{idx}_thickness"] = float(step.get('thickness', 10.0)) if 'bias_slider_val' in config: st.session_state['bias_slider_val'] = float(config['bias_slider_val']) if 'enable_avalanche_toggle_val' in config: st.session_state['enable_avalanche_toggle_val'] = bool(config['enable_avalanche_toggle_val']) + if 'enable_btbt_toggle_val' in config: + st.session_state['enable_btbt_toggle_val'] = bool(config['enable_btbt_toggle_val']) if 'run_full_sweep_toggle_val' in config: st.session_state['run_full_sweep_toggle_val'] = bool(config['run_full_sweep_toggle_val']) if 'plot_doping_xmax' in config: st.session_state['plot_doping_xmax'] = float(config['plot_doping_xmax']) if 'plot_electro_xmax' in config: st.session_state['plot_electro_xmax'] = float(config['plot_electro_xmax']) + if 'sweep_v_max' in config: + st.session_state['sweep_v_max'] = float(config['sweep_v_max']) st.session_state['last_loaded_config_hash'] = file_hash st.success("Config retrieved successfully!") @@ -395,15 +403,25 @@ with st.sidebar: min_value=0.1, max_value=500.0, value=5.0, step=0.5, help="Electrostatic Profile X-axis maximum range (μm)." ) + sweep_v_max = label_input_row( + "Sweep Max (V)", "number_input", "sweep_v_max", + ratio=[7, 3], + min_value=1.0, max_value=1000.0, value=1000.0, step=10.0, + help="Sweep target voltage (e.g. 1000V for high-voltage, 15V for Zener)." + ) st.markdown("### ⚡ Bias & Physics") + sweep_max_v = float(st.session_state.get('sweep_v_max', 1000.0)) + if st.session_state.get('bias_slider_val', 5.0) > sweep_max_v: + st.session_state['bias_slider_val'] = sweep_max_v + bias = label_input_row( "Bias Vbias (V)", "slider", "bias_slider_val", ratio=[4.5, 5.5], min_value=0.0, - max_value=1000.0, - value=5.0, - step=1.0, + max_value=sweep_max_v, + value=min(5.0, sweep_max_v), + step=0.1 if sweep_max_v <= 20.0 else 1.0, help="Applied at bottom contact (x = L). Top contact (x = 0) is reference 0V." ) @@ -414,11 +432,18 @@ with st.sidebar: help="Enable impact ionization for the selected voltage simulation." ) + enable_btbt = label_input_row( + "BTBT (Tunneling)", "toggle", "enable_btbt_toggle_val", + ratio=[6.5, 3.5], + value=False, + help="Enable band-to-band tunneling for the selected voltage simulation." + ) + run_full_sweep = label_input_row( - "Run 1000V Sweep", "toggle", "run_full_sweep_toggle_val", + f"Run {sweep_max_v:.0f}V Sweep", "toggle", "run_full_sweep_toggle_val", ratio=[6.5, 3.5], value=True, - help="Run the full 0~1000V sweep with and without avalanche. Turn off for faster parameter tuning." + help=f"Run the full 0~{sweep_max_v:.0f}V sweep with and without avalanche. Turn off for faster parameter tuning." ) st.markdown("---") @@ -464,6 +489,7 @@ with st.sidebar: 'energy': 80.0, 'dose': 1e12, 'cs': 1e19, + 'thickness': 10.0, 'temp': 1000.0, 'time': 60.0 } @@ -481,6 +507,7 @@ with st.sidebar: 'energy': 80.0, 'dose': 1e12, 'cs': 1e19, + 'thickness': 10.0, 'temp': 1000.0, 'time': 60.0 }) @@ -539,7 +566,7 @@ with st.sidebar: dopant_idx = dopant_options.index(step['dopant']) if step['dopant'] in dopant_options else 0 dopant = label_input_row("Dopant", "selectbox", f"step_{idx}_dopant", options=dopant_options, index=dopant_idx) - method_options = ['Implant', 'Diffusion'] + method_options = ['Implant', 'Diffusion', 'Epi Growth'] method_idx = method_options.index(step['method']) if step['method'] in method_options else 0 method = label_input_row("Method", "selectbox", f"step_{idx}_method", options=method_options, index=method_idx) @@ -548,13 +575,21 @@ with st.sidebar: energy = label_input_row("Energy (keV)", "number_input", f"step_{idx}_energy", ratio=[7, 3], min_value=1.0, max_value=1000.0, value=float(step.get('energy', 80.0)), step=5.0) dose = label_input_row("Dose (cm⁻²)", "number_input", f"step_{idx}_dose", ratio=[7, 3], min_value=1e9, max_value=1e17, value=float(step.get('dose', 1e12)), format="%e") cs = step.get('cs', 1e19) # preserve - else: # Diffusion + thickness = step.get('thickness', 10.0) # preserve + elif method == 'Diffusion': # Diffusion surface_options = ['Top', 'Bottom'] surface_idx = surface_options.index(step.get('surface', 'Top')) if step.get('surface', 'Top') in surface_options else 0 surface = label_input_row("Surface Loc.", "selectbox", f"step_{idx}_surface", ratio=[5, 5], options=surface_options, index=surface_idx) cs = label_input_row("Surface Cs (cm⁻³)", "number_input", f"step_{idx}_cs", ratio=[7, 3], min_value=1e13, max_value=1e22, value=float(step.get('cs', 1e19)), format="%e") energy = step.get('energy', 80.0) # preserve dose = step.get('dose', 1e12) # preserve + thickness = step.get('thickness', 10.0) # preserve + else: # Epi Growth + surface = 'Top' + thickness = label_input_row("Thickness (μm)", "number_input", f"step_{idx}_thickness", ratio=[7, 3], min_value=0.1, max_value=200.0, value=float(step.get('thickness', 10.0)), step=0.5) + cs = label_input_row("Doping Conc. (cm⁻³)", "number_input", f"step_{idx}_cs", ratio=[7, 3], min_value=1e10, max_value=1e20, value=float(step.get('cs', 1e15)), format="%e") + energy = step.get('energy', 80.0) # preserve + dose = step.get('dose', 1e12) # preserve temp = label_input_row("Temp (°C)", "number_input", f"step_{idx}_temp", ratio=[7, 3], min_value=25.0, max_value=1300.0, value=float(step['temp']), step=25.0) time = label_input_row("Time (min)", "number_input", f"step_{idx}_time", ratio=[7, 3], min_value=0.0, max_value=1000.0, value=float(step['time']), step=5.0) @@ -568,6 +603,7 @@ with st.sidebar: 'energy': energy, 'dose': dose, 'cs': cs, + 'thickness': thickness, 'temp': temp, 'time': time }) @@ -600,6 +636,7 @@ with st.sidebar: 'energy': 80.0, 'dose': 1e12, 'cs': 1e19, + 'thickness': 10.0, 'temp': 1000.0, 'time': 60.0 } @@ -617,9 +654,11 @@ with st.sidebar: 'process_steps': st.session_state.get('process_steps', []), 'bias_slider_val': st.session_state.get('bias_slider_val', 5.0), 'enable_avalanche_toggle_val': st.session_state.get('enable_avalanche_toggle_val', False), + 'enable_btbt_toggle_val': st.session_state.get('enable_btbt_toggle_val', False), 'run_full_sweep_toggle_val': st.session_state.get('run_full_sweep_toggle_val', True), 'plot_doping_xmax': st.session_state.get('plot_doping_xmax', 5.0), 'plot_electro_xmax': st.session_state.get('plot_electro_xmax', 5.0), + 'sweep_v_max': st.session_state.get('sweep_v_max', 1000.0), } json_str = json.dumps(config_to_save, indent=2) download_btn_container.download_button( @@ -643,27 +682,48 @@ st.markdown( # --- 6. Execute Simulation & Cache --- # Create unique key to track doping state changes -doping_key = (substrate_type, n_sub, length, area_cm2, str(st.session_state['process_steps'])) +sweep_max_v = float(st.session_state.get('sweep_v_max', 1000.0)) +doping_key = (substrate_type, n_sub, length, area_cm2, str(st.session_state['process_steps']), sweep_max_v) # Initialize I-V curve caches if not present +if 'iv_curve_basic' not in st.session_state: + st.session_state['iv_curve_basic'] = ([0.0], [0.0]) if 'iv_curve_with_avalanche' not in st.session_state: st.session_state['iv_curve_with_avalanche'] = ([0.0], [0.0]) -if 'iv_curve_without_avalanche' not in st.session_state: - st.session_state['iv_curve_without_avalanche'] = ([0.0], [0.0]) +if 'iv_curve_with_btbt' not in st.session_state: + st.session_state['iv_curve_with_btbt'] = ([0.0], [0.0]) try: with DEVSIM_LOCK: if run_full_sweep and ('cached_doping_key' not in st.session_state or st.session_state['cached_doping_key'] != doping_key): - with st.spinner("Calculating high-voltage I-V sweeps (0 ~ 1000V)..."): - # 1. Sweep with avalanche + with st.spinner(f"Calculating I-V sweeps (0 ~ {sweep_max_v:.0f}V)..."): + # 1. Sweep Basic (no avalanche, no btbt) + try: + res_no_av = build_and_solve_1d( + bias_target=sweep_max_v, + substrate_type=substrate_type, + substrate_doping=n_sub, + length=length, + process_steps=st.session_state['process_steps'], + enable_avalanche=False, + enable_btbt=False, + area_cm2=area_cm2 + ) + v_no_av, j_no_av = res_no_av['v_history'], res_no_av['j_history'] + except Exception as e: + v_no_av, j_no_av = [0.0], [0.0] + st.error(f"Basic sweep failed: {e}") + + # 2. Sweep with avalanche (no btbt) try: res_av = build_and_solve_1d( - bias_target=1000.0, + bias_target=sweep_max_v, substrate_type=substrate_type, substrate_doping=n_sub, length=length, process_steps=st.session_state['process_steps'], enable_avalanche=True, + enable_btbt=False, area_cm2=area_cm2 ) v_av, j_av = res_av['v_history'], res_av['j_history'] @@ -671,28 +731,38 @@ try: v_av, j_av = [0.0], [0.0] st.error(f"Avalanche sweep failed to converge: {e}") - # 2. Sweep without avalanche + # 3. Sweep with BTBT (no avalanche) try: - res_no_av = build_and_solve_1d( - bias_target=1000.0, + res_btbt = build_and_solve_1d( + bias_target=sweep_max_v, substrate_type=substrate_type, substrate_doping=n_sub, length=length, process_steps=st.session_state['process_steps'], enable_avalanche=False, + enable_btbt=True, area_cm2=area_cm2 ) - v_no_av, j_no_av = res_no_av['v_history'], res_no_av['j_history'] + v_btbt, j_btbt = res_btbt['v_history'], res_btbt['j_history'] except Exception as e: - v_no_av, j_no_av = [0.0], [0.0] - st.error(f"Without-avalanche sweep failed: {e}") + v_btbt, j_btbt = [0.0], [0.0] + st.error(f"BTBT sweep failed to converge: {e}") # Store in session state st.session_state['cached_doping_key'] = doping_key + st.session_state['iv_curve_basic'] = (v_no_av, j_no_av) st.session_state['iv_curve_with_avalanche'] = (v_av, j_av) - st.session_state['iv_curve_without_avalanche'] = (v_no_av, j_no_av) + st.session_state['iv_curve_with_btbt'] = (v_btbt, j_btbt) + v_pt_val = None + if 'res_no_av' in locals() and isinstance(res_no_av, dict): + v_pt_val = res_no_av.get('v_punchthrough') + if v_pt_val is None and 'res_av' in locals() and isinstance(res_av, dict): + v_pt_val = res_av.get('v_punchthrough') + if v_pt_val is None and 'res_btbt' in locals() and isinstance(res_btbt, dict): + v_pt_val = res_btbt.get('v_punchthrough') + st.session_state['v_punchthrough'] = v_pt_val - # 3. Single-bias simulation for electrostatic plots + # 4. Single-bias simulation for electrostatic plots with st.spinner(f"Solving electrostatic profiles at Vbias = {bias} V..."): res = build_and_solve_1d( bias_target=bias, @@ -701,8 +771,12 @@ try: length=length, process_steps=st.session_state['process_steps'], enable_avalanche=enable_avalanche, + enable_btbt=enable_btbt, area_cm2=area_cm2 ) + # If full sweep did not run or v_pt was not found, try to grab from single bias solve + if not run_full_sweep or st.session_state.get('v_punchthrough') is None: + st.session_state['v_punchthrough'] = res.get('v_punchthrough') # Check if target bias was reached v_actual = res["v_solved"] @@ -850,7 +924,23 @@ try: with col_right: figC = go.Figure() - # 1. Curve with avalanche + # 1. Curve: Basic (Without Avalanche, Without BTBT) + v_basic, j_basic = st.session_state.get('iv_curve_basic', ([0.0], [0.0])) + if len(v_basic) == 1 and v_basic[0] == 0.0: + # Fallback to legacy cached key + v_basic, j_basic = st.session_state.get('iv_curve_without_avalanche', ([0.0], [0.0])) + i_basic_ma = np.abs(np.array(j_basic)) * area_cm2 * 1000.0 + i_basic_ma = np.clip(i_basic_ma, 1e-12, None) + figC.add_trace(go.Scatter( + x=v_basic, y=i_basic_ma, + name="Basic (Drift-Diffusion)", + line=dict(color='#00d2ff', width=2, dash='dot'), + mode='lines+markers', + marker=dict(size=4), + hovertemplate='Bias: %{x:.1f} V
Current: %{y:.3e} mA' + )) + + # 2. Curve: With Avalanche v_av, j_av = st.session_state.get('iv_curve_with_avalanche', ([0.0], [0.0])) i_av_ma = np.abs(np.array(j_av)) * area_cm2 * 1000.0 i_av_ma = np.clip(i_av_ma, 1e-12, None) @@ -862,21 +952,34 @@ try: hovertemplate='Bias: %{x:.1f} V
Current: %{y:.3e} mA' )) - # 2. Curve without avalanche - v_no_av, j_no_av = st.session_state.get('iv_curve_without_avalanche', ([0.0], [0.0])) - i_no_av_ma = np.abs(np.array(j_no_av)) * area_cm2 * 1000.0 - i_no_av_ma = np.clip(i_no_av_ma, 1e-12, None) + # 3. Curve: With BTBT + v_btbt, j_btbt = st.session_state.get('iv_curve_with_btbt', ([0.0], [0.0])) + i_btbt_ma = np.abs(np.array(j_btbt)) * area_cm2 * 1000.0 + i_btbt_ma = np.clip(i_btbt_ma, 1e-12, None) figC.add_trace(go.Scatter( - x=v_no_av, y=i_no_av_ma, - name="Without Avalanche", - line=dict(color='#00d2ff', width=2, dash='dot'), - mode='lines', + x=v_btbt, y=i_btbt_ma, + name="With BTBT", + line=dict(color='#a000ff', width=3), + mode='lines+markers', hovertemplate='Bias: %{x:.1f} V
Current: %{y:.3e} mA' )) + # Add punch-through vertical line if available + v_pt = st.session_state.get('v_punchthrough') + if v_pt is not None: + figC.add_vline( + x=v_pt, + line_width=2, + line_dash="dash", + line_color="#ff4b4b", + annotation_text=f"Punch-through: {v_pt:.1f}V", + annotation_position="top left", + annotation_font=dict(color="#ff4b4b", size=12) + ) + figC.update_layout(**plot_layout) figC.update_layout( - title="📉 High-Voltage I-V Characteristics (0 ~ 1000V)", + title=f"📉 I-V Characteristics (0 ~ {sweep_max_v:.0f}V)", height=750, legend=dict(y=-0.12), xaxis=dict(title="Applied Bias Vbias (V)"), diff --git a/gui1d/sim1d_config (2).json b/gui1d/sim1d_config (2).json new file mode 100644 index 0000000..7539453 --- /dev/null +++ b/gui1d/sim1d_config (2).json @@ -0,0 +1,37 @@ +{ + "sub_type": "n", + "sub_doping": 55000000000000.0, + "sub_length": 100.0, + "device_area": 0.01, + "process_steps": [ + { + "enabled": true, + "type": "p", + "dopant": "Boron", + "method": "Implant", + "surface": "Top", + "energy": 35.0, + "dose": 6000000000000.0, + "cs": 1e+19, + "temp": 1000.0, + "time": 60.0 + }, + { + "enabled": true, + "type": "p", + "dopant": "Boron", + "method": "Implant", + "surface": "Top", + "energy": 90.0, + "dose": 300000000000000.0, + "cs": 1e+19, + "temp": 1150.0, + "time": 360.0 + } + ], + "bias_slider_val": 188.0, + "enable_avalanche_toggle_val": false, + "run_full_sweep_toggle_val": true, + "plot_doping_xmax": 16.0, + "plot_electro_xmax": 100.0 +} \ No newline at end of file diff --git a/gui1d/sim1d_config (3).json b/gui1d/sim1d_config (3).json new file mode 100644 index 0000000..8a212b9 --- /dev/null +++ b/gui1d/sim1d_config (3).json @@ -0,0 +1,37 @@ +{ + "sub_type": "n", + "sub_doping": 55000000000000.0, + "sub_length": 200.0, + "device_area": 0.01, + "process_steps": [ + { + "enabled": true, + "type": "p", + "dopant": "Boron", + "method": "Implant", + "surface": "Top", + "energy": 35.0, + "dose": 6000000000000.0, + "cs": 1e+19, + "temp": 1000.0, + "time": 60.0 + }, + { + "enabled": true, + "type": "p", + "dopant": "Boron", + "method": "Implant", + "surface": "Top", + "energy": 90.0, + "dose": 300000000000000.0, + "cs": 1e+19, + "temp": 1150.0, + "time": 360.0 + } + ], + "bias_slider_val": 953.0, + "enable_avalanche_toggle_val": false, + "run_full_sweep_toggle_val": true, + "plot_doping_xmax": 16.0, + "plot_electro_xmax": 200.0 +} \ No newline at end of file diff --git a/gui1d/sim1d_config5Vzener.json b/gui1d/sim1d_config5Vzener.json new file mode 100644 index 0000000..588f6ba --- /dev/null +++ b/gui1d/sim1d_config5Vzener.json @@ -0,0 +1,26 @@ +{ + "sub_type": "n", + "sub_doping": 2e+18, + "sub_length": 100.0, + "device_area": 0.01, + "process_steps": [ + { + "enabled": true, + "type": "p", + "dopant": "Boron", + "method": "Diffusion", + "surface": "Top", + "energy": 80.0, + "dose": 1000000000000.0, + "cs": 1e+19, + "temp": 1000.0, + "time": 30.0 + } + ], + "bias_slider_val": 5.0, + "enable_avalanche_toggle_val": true, + "enable_btbt_toggle_val": true, + "run_full_sweep_toggle_val": true, + "plot_doping_xmax": 1.5, + "plot_electro_xmax": 1.5 +} \ No newline at end of file diff --git a/gui1d/sim1d_config5Vzener2.json b/gui1d/sim1d_config5Vzener2.json new file mode 100644 index 0000000..9320cd2 --- /dev/null +++ b/gui1d/sim1d_config5Vzener2.json @@ -0,0 +1,27 @@ +{ + "sub_type": "n", + "sub_doping": 2e+18, + "sub_length": 10.0, + "device_area": 0.01, + "process_steps": [ + { + "enabled": true, + "type": "p", + "dopant": "Boron", + "method": "Diffusion", + "surface": "Top", + "energy": 80.0, + "dose": 1000000000000.0, + "cs": 1e+19, + "temp": 1000.0, + "time": 15.0 + } + ], + "bias_slider_val": 5.0, + "enable_avalanche_toggle_val": true, + "enable_btbt_toggle_val": true, + "run_full_sweep_toggle_val": true, + "plot_doping_xmax": 1.5, + "plot_electro_xmax": 1.5, + "sweep_v_max": 10.0 +} \ No newline at end of file diff --git a/gui1d/sim1d_config800V.json b/gui1d/sim1d_config800V.json new file mode 100644 index 0000000..5f881c2 --- /dev/null +++ b/gui1d/sim1d_config800V.json @@ -0,0 +1,37 @@ +{ + "sub_type": "n", + "sub_doping": 200000000000000.0, + "sub_length": 200.0, + "device_area": 0.01, + "process_steps": [ + { + "enabled": true, + "type": "p", + "dopant": "Boron", + "method": "Implant", + "surface": "Top", + "energy": 35.0, + "dose": 6000000000000.0, + "cs": 1e+19, + "temp": 1000.0, + "time": 60.0 + }, + { + "enabled": true, + "type": "p", + "dopant": "Boron", + "method": "Implant", + "surface": "Top", + "energy": 90.0, + "dose": 300000000000000.0, + "cs": 1e+19, + "temp": 1150.0, + "time": 360.0 + } + ], + "bias_slider_val": 953.0, + "enable_avalanche_toggle_val": false, + "run_full_sweep_toggle_val": true, + "plot_doping_xmax": 16.0, + "plot_electro_xmax": 200.0 +} \ No newline at end of file diff --git a/gui1d/sim1d_config_NBL1.json b/gui1d/sim1d_config_NBL1.json new file mode 100644 index 0000000..d772f76 --- /dev/null +++ b/gui1d/sim1d_config_NBL1.json @@ -0,0 +1,54 @@ +{ + "sub_type": "n", + "sub_doping": 1000000000000000.0, + "sub_length": 50.0, + "device_area": 0.01, + "process_steps": [ + { + "enabled": true, + "type": "n", + "dopant": "Phosphorus", + "method": "Implant", + "surface": "Top", + "energy": 80.0, + "dose": 5000000000000000.0, + "cs": 1000000000000000.0, + "thickness": 10.0, + "temp": 1150.0, + "time": 180.0 + }, + { + "enabled": true, + "type": "n", + "dopant": "Arsenic", + "method": "Epi Growth", + "surface": "Top", + "energy": 80.0, + "dose": 1000000000000.0, + "cs": 1000000000000000.0, + "thickness": 10.0, + "temp": 1050.0, + "time": 30.0 + }, + { + "enabled": true, + "type": "p", + "dopant": "Boron", + "method": "Implant", + "surface": "Top", + "energy": 30.0, + "dose": 1000000000000.0, + "cs": 1e+19, + "thickness": 10.0, + "temp": 950.0, + "time": 30.0 + } + ], + "bias_slider_val": 56.0, + "enable_avalanche_toggle_val": false, + "enable_btbt_toggle_val": false, + "run_full_sweep_toggle_val": true, + "plot_doping_xmax": 30.0, + "plot_electro_xmax": 30.0, + "sweep_v_max": 1000.0 +} \ No newline at end of file diff --git a/gui1d/sim1d_config_NBL1a.json b/gui1d/sim1d_config_NBL1a.json new file mode 100644 index 0000000..fa889f1 --- /dev/null +++ b/gui1d/sim1d_config_NBL1a.json @@ -0,0 +1,54 @@ +{ + "sub_type": "n", + "sub_doping": 1000000000000000.0, + "sub_length": 50.0, + "device_area": 0.01, + "process_steps": [ + { + "enabled": true, + "type": "n", + "dopant": "Phosphorus", + "method": "Implant", + "surface": "Top", + "energy": 80.0, + "dose": 5000000000000000.0, + "cs": 1000000000000000.0, + "thickness": 10.0, + "temp": 1150.0, + "time": 60.0 + }, + { + "enabled": true, + "type": "n", + "dopant": "Arsenic", + "method": "Epi Growth", + "surface": "Top", + "energy": 80.0, + "dose": 1000000000000.0, + "cs": 1000000000000000.0, + "thickness": 10.0, + "temp": 1050.0, + "time": 30.0 + }, + { + "enabled": true, + "type": "p", + "dopant": "Boron", + "method": "Implant", + "surface": "Top", + "energy": 30.0, + "dose": 1000000000000.0, + "cs": 1e+19, + "thickness": 10.0, + "temp": 950.0, + "time": 30.0 + } + ], + "bias_slider_val": 30.0, + "enable_avalanche_toggle_val": false, + "enable_btbt_toggle_val": false, + "run_full_sweep_toggle_val": true, + "plot_doping_xmax": 30.0, + "plot_electro_xmax": 30.0, + "sweep_v_max": 1000.0 +} \ No newline at end of file diff --git a/gui1d/sim1d_config_bd&pt.json b/gui1d/sim1d_config_bd&pt.json new file mode 100644 index 0000000..8bd727f --- /dev/null +++ b/gui1d/sim1d_config_bd&pt.json @@ -0,0 +1,26 @@ +{ + "sub_type": "n", + "sub_doping": 3e+16, + "sub_length": 200.0, + "device_area": 0.01, + "process_steps": [ + { + "enabled": true, + "type": "p", + "dopant": "Boron", + "method": "Implant", + "surface": "Top", + "energy": 35.0, + "dose": 6000000000000.0, + "cs": 1e+19, + "temp": 1000.0, + "time": 60.0 + } + ], + "bias_slider_val": 953.0, + "enable_avalanche_toggle_val": true, + "enable_btbt_toggle_val": true, + "run_full_sweep_toggle_val": true, + "plot_doping_xmax": 5.0, + "plot_electro_xmax": 5.0 +} \ No newline at end of file diff --git a/gui1d/sim1d_config_bd&pt1.json b/gui1d/sim1d_config_bd&pt1.json new file mode 100644 index 0000000..7d31956 --- /dev/null +++ b/gui1d/sim1d_config_bd&pt1.json @@ -0,0 +1,26 @@ +{ + "sub_type": "n", + "sub_doping": 1e+17, + "sub_length": 200.0, + "device_area": 0.01, + "process_steps": [ + { + "enabled": true, + "type": "p", + "dopant": "Boron", + "method": "Diffusion", + "surface": "Top", + "energy": 35.0, + "dose": 6000000000000.0, + "cs": 1e+19, + "temp": 1100.0, + "time": 180.0 + } + ], + "bias_slider_val": 19.0, + "enable_avalanche_toggle_val": true, + "enable_btbt_toggle_val": true, + "run_full_sweep_toggle_val": true, + "plot_doping_xmax": 5.0, + "plot_electro_xmax": 5.0 +} \ No newline at end of file diff --git a/gui1d/sim1d_config_punch725V.json b/gui1d/sim1d_config_punch725V.json new file mode 100644 index 0000000..2505b1f --- /dev/null +++ b/gui1d/sim1d_config_punch725V.json @@ -0,0 +1,25 @@ +{ + "sub_type": "n", + "sub_doping": 55000000000000.0, + "sub_length": 200.0, + "device_area": 0.01, + "process_steps": [ + { + "enabled": true, + "type": "p", + "dopant": "Boron", + "method": "Implant", + "surface": "Top", + "energy": 35.0, + "dose": 700000000000.0, + "cs": 1e+19, + "temp": 1000.0, + "time": 60.0 + } + ], + "bias_slider_val": 953.0, + "enable_avalanche_toggle_val": true, + "run_full_sweep_toggle_val": true, + "plot_doping_xmax": 16.0, + "plot_electro_xmax": 200.0 +} \ No newline at end of file diff --git a/gui1d/solve_1d.py b/gui1d/solve_1d.py index 2d5cf28..988dcbc 100644 --- a/gui1d/solve_1d.py +++ b/gui1d/solve_1d.py @@ -18,7 +18,8 @@ from physics.new_physics import ( CreateHFMobility, CreateSiliconDriftDiffusion, CreateSiliconDriftDiffusionContact, - CreateAvalancheGeneration + CreateAvalancheGeneration, + CreateBTBTGeneration ) # Vectorized complementary error function @@ -74,16 +75,22 @@ def calc_donors_acceptors(x_um, process_steps, substrate_type, substrate_doping, """ Calculate the separate Donors and Acceptors concentration arrays (cm^-3) across a position array x_um (in microns), and also return individual step profiles. + length here is total_length. """ donors = np.zeros_like(x_um) acceptors = np.zeros_like(x_um) step_profiles = [] + # Calculate total epi thickness to know where substrate is + epi_thickness = sum(step.get('thickness', 0.0) for step in process_steps if step.get('enabled', True) and step.get('method') == 'Epi Growth') + substrate_thickness = length - epi_thickness + # Initialize substrate baseline + # Substrate occupies [epi_thickness, length] if substrate_type == 'n': - donors += substrate_doping + donors += np.where(x_um >= epi_thickness, substrate_doping, 0.0) else: - acceptors += substrate_doping + acceptors += np.where(x_um >= epi_thickness, substrate_doping, 0.0) # Add background contact doping to ensure Ohmic contact at bottom (x = L) # Contact doping thickness is 0.5 um, peak is 1e19 @@ -97,7 +104,7 @@ def calc_donors_acceptors(x_um, process_steps, substrate_type, substrate_doping, step_profiles.append({ 'name': 'Substrate (Step 0)', 'type': substrate_type, - 'profile': np.full_like(x_um, substrate_doping) + 'profile': np.where(x_um >= epi_thickness, substrate_doping, 0.0) }) # Compute profiles for each process step @@ -126,6 +133,10 @@ def calc_donors_acceptors(x_um, process_steps, substrate_type, substrate_doping, # Convert Dt to microns^2 Dt_total_um2 = Dt_total_cm2 * 1e8 + # Calculate reference position for the surface at step i. + # It is shifted by all Epi Growth steps that occur AFTER step i. + x_left = sum(sub_step.get('thickness', 0.0) for sub_step in process_steps[i+1:] if sub_step.get('enabled', True) and sub_step.get('method') == 'Epi Growth') + # Calculate profile if method == 'Implant': energy = step.get('energy', 80.0) @@ -135,22 +146,34 @@ def calc_donors_acceptors(x_um, process_steps, substrate_type, substrate_doping, peak_conc = dose / (np.sqrt(2 * np.pi) * (dRp_eff * 1e-4)) # dose in cm^-2, dRp in cm if surface == 'Top': - prof = peak_conc * (np.exp(-((x_um - Rp) / (np.sqrt(2) * dRp_eff)) ** 2) + - np.exp(-((x_um + Rp) / (np.sqrt(2) * dRp_eff)) ** 2)) + prof = peak_conc * (np.exp(-((x_um - x_left - Rp) / (np.sqrt(2) * dRp_eff)) ** 2) + + np.exp(-((x_um - x_left + Rp) / (np.sqrt(2) * dRp_eff)) ** 2)) else: prof = peak_conc * (np.exp(-((x_um - (length - Rp)) / (np.sqrt(2) * dRp_eff)) ** 2) + np.exp(-((x_um - (length + Rp)) / (np.sqrt(2) * dRp_eff)) ** 2)) - else: # Constant Source Predeposition + elif method == 'Diffusion': # Constant Source Predeposition cs = step.get('cs', 1e19) if Dt_total_um2 <= 0.0: Dt_total_um2 = 1e-10 if surface == 'Top': - prof = cs * erfc_vec(x_um / (2 * np.sqrt(Dt_total_um2))) + prof = cs * erfc_vec((x_um - x_left) / (2 * np.sqrt(Dt_total_um2))) else: prof = cs * erfc_vec((length - x_um) / (2 * np.sqrt(Dt_total_um2))) + else: # Epi Growth + thick_val = step.get('thickness', 10.0) + cs_val = step.get('cs', 1e15) + x_right = x_left + thick_val + + if Dt_total_um2 <= 0.0: + Dt_total_um2 = 1e-10 + + # Slab diffusion formula + prof = 0.5 * cs_val * (erfc_vec((x_um - x_right) / (2 * np.sqrt(Dt_total_um2))) - + erfc_vec((x_um - x_left) / (2 * np.sqrt(Dt_total_um2)))) + # Add to Net Donors/Acceptors arrays if type_ == 'n': donors += prof @@ -183,7 +206,6 @@ def find_junction_depths(process_steps, substrate_type, substrate_doping, length junctions.append(x_cross) return junctions - def build_and_solve_1d( bias_target, substrate_type='n', @@ -191,12 +213,17 @@ def build_and_solve_1d( length=30.0, process_steps=[], enable_avalanche=False, + enable_btbt=False, area_cm2=1.0 ): """ Builds a 1D Diode mesh, sets up doping, solves from 0V equilibrium to bias_target using a rapid micro-sweep, and returns physical profiles. """ + # Calculate total thickness including Epi Growth steps + epi_thickness = sum(step.get('thickness', 0.0) for step in process_steps if step.get('enabled', True) and step.get('method') == 'Epi Growth') + total_length = length + epi_thickness + # 1. Reset DEVSIM state to prevent name collisions devsim.reset_devsim() @@ -208,13 +235,20 @@ def build_and_solve_1d( # 2. Build 1D adaptive mesh # Find junctions to refine the mesh around them - junctions = find_junction_depths(process_steps, substrate_type, substrate_doping, length) + junctions = find_junction_depths(process_steps, substrate_type, substrate_doping, total_length) # Define control points (x, target spacing) + # Higher doping -> narrower depletion region -> needs finer mesh (down to 0.5nm) + junction_spacing = max(0.0005, min(0.02, 4e15 / float(substrate_doping))) + control_points = [(0.0, 0.05)] for j in junctions: - control_points.append((j, 0.02)) # Fine mesh around junctions - control_points.append((length, 0.5)) # Coarser mesh near bottom + # Refine a symmetric 100nm region around the junction to ensure + # smooth electric field profiles even when depletion region expands under bias + control_points.append((max(0.0, j - 0.05), junction_spacing)) + control_points.append((j, junction_spacing)) + control_points.append((min(total_length, j + 0.05), junction_spacing)) + control_points.append((total_length, 0.5)) # Coarser mesh near bottom # Sort control points control_points.sort(key=lambda x: x[0]) @@ -232,7 +266,7 @@ def build_and_solve_1d( while curr_x < x_end: t = (curr_x - x_start) / seg_len spacing = sp_start + t * (sp_end - sp_start) - spacing = max(0.01, min(2.0, spacing)) + spacing = max(0.0005, min(2.0, spacing)) curr_x += spacing if curr_x < x_end - 0.005: all_x.append(curr_x) @@ -262,12 +296,15 @@ def build_and_solve_1d( devsim.add_1d_contact(mesh=mesh_name, name="bottom", tag="bottom", material="metal") devsim.add_1d_region(mesh=mesh_name, region=region, tag1="top", tag2="bottom", material="Silicon") devsim.finalize_mesh(mesh=mesh_name) - + devsim.create_device(mesh=mesh_name, device=device) # 3. Setup Doping Profile models in DEVSIM - # Calculate separate donors and acceptors on the grid - donors_array, acceptors_array, step_profiles = calc_donors_acceptors(all_x, process_steps, substrate_type, substrate_doping, length) + # Get actual finalized node coordinates from DEVSIM (guarantees length matching) + x_devsim = np.array(devsim.get_node_model_values(device=device, region=region, name="x")) / um + + # Calculate separate donors and acceptors on the actual grid coordinates + donors_array, acceptors_array, step_profiles = calc_donors_acceptors(x_devsim, process_steps, substrate_type, substrate_doping, total_length) # Set Donors directly devsim.node_solution(device=device, region=region, name="Donors_data") @@ -321,25 +358,67 @@ def build_and_solve_1d( # Override equations if enable_avalanche: CreateAvalancheGeneration(device, region, hf_opts['Jn'], hf_opts['Jp']) - devsim.equation(device=device, region=region, name="ElectronContinuityEquation", variable_name="Electrons", - time_node_model="NCharge", edge_model=hf_opts['Jn'], edge_volume_model="AvalancheGeneration", - variable_update="positive", node_model="ElectronGeneration", min_error=1e5) - devsim.equation(device=device, region=region, name="HoleContinuityEquation", variable_name="Holes", - time_node_model="PCharge", edge_model=hf_opts['Jp'], edge_volume_model="AvalancheGeneration_p", - variable_update="positive", node_model="HoleGeneration", min_error=1e5) + if enable_btbt: + CreateBTBTGeneration(device, region) + + av_model_n = "AvalancheGeneration" if enable_avalanche else "" + av_model_p = "AvalancheGeneration_p" if enable_avalanche else "" + btbt_model_n = "BTBTGeneration" if enable_btbt else "" + btbt_model_p = "BTBTGeneration_p" if enable_btbt else "" + + if av_model_n and btbt_model_n: + from physics.model_create import CreateEdgeModel, CreateEdgeModelDerivatives + CreateEdgeModel(device, region, "CombinedGeneration", "AvalancheGeneration + BTBTGeneration") + CreateEdgeModel(device, region, "CombinedGeneration_p", "AvalancheGeneration_p + BTBTGeneration_p") + for i in ("Potential", "Electrons", "Holes"): + CreateEdgeModelDerivatives(device, region, "CombinedGeneration", "AvalancheGeneration + BTBTGeneration", i) + CreateEdgeModelDerivatives(device, region, "CombinedGeneration_p", "AvalancheGeneration_p + BTBTGeneration_p", i) + gen_model_n = "CombinedGeneration" + gen_model_p = "CombinedGeneration_p" + elif av_model_n: + gen_model_n = av_model_n + gen_model_p = av_model_p + elif btbt_model_n: + gen_model_n = btbt_model_n + gen_model_p = btbt_model_p else: - devsim.equation(device=device, region=region, name="ElectronContinuityEquation", variable_name="Electrons", - time_node_model="NCharge", edge_model=hf_opts['Jn'], - variable_update="positive", node_model="ElectronGeneration", min_error=1e5) - devsim.equation(device=device, region=region, name="HoleContinuityEquation", variable_name="Holes", - time_node_model="PCharge", edge_model=hf_opts['Jp'], - variable_update="positive", node_model="HoleGeneration", min_error=1e5) - + gen_model_n = "" + gen_model_p = "" + + # First set up equations WITHOUT generation terms (basic drift-diffusion) + devsim.equation(device=device, region=region, name="ElectronContinuityEquation", variable_name="Electrons", + time_node_model="NCharge", edge_model=hf_opts['Jn'], + variable_update="positive", node_model="ElectronGeneration", min_error=1e5) + devsim.equation(device=device, region=region, name="HoleContinuityEquation", variable_name="Holes", + time_node_model="PCharge", edge_model=hf_opts['Jp'], + variable_update="positive", node_model="HoleGeneration", min_error=1e5) devsim.equation(device=device, region=region, name="PotentialEquation", variable_name="Potential", node_model="PotentialNodeCharge", edge_model="DField", variable_update="default", min_error=1e-3) - # Solve 0V Drift-Diffusion - devsim.solve(type="dc", absolute_error=1e10, relative_error=1e-8, charge_error=1e12, maximum_iterations=100) + # Solve 0V basic Drift-Diffusion (highly stable) + devsim.solve(type="dc", absolute_error=1e10, relative_error=1e-5, charge_error=1e12, maximum_iterations=100) + + # Now if generation models are enabled, add them and re-solve at 0V starting from the solved state + if gen_model_n: + devsim.equation(device=device, region=region, name="ElectronContinuityEquation", variable_name="Electrons", + time_node_model="NCharge", edge_model=hf_opts['Jn'], edge_volume_model=gen_model_n, + variable_update="positive", node_model="ElectronGeneration", min_error=1e5) + devsim.equation(device=device, region=region, name="HoleContinuityEquation", variable_name="Holes", + time_node_model="PCharge", edge_model=hf_opts['Jp'], edge_volume_model=gen_model_p, + variable_update="positive", node_model="HoleGeneration", min_error=1e5) + # Re-solve at 0V to incorporate generation models smoothly + devsim.solve(type="dc", absolute_error=1e10, relative_error=1e-5, charge_error=1e12, maximum_iterations=100) + + # Check 0V depletion for punch-through + electrons_0 = np.array(devsim.get_node_model_values(device=device, region=region, name="Electrons")) + holes_0 = np.array(devsim.get_node_model_values(device=device, region=region, name="Holes")) + doping_0 = np.array(devsim.get_node_model_values(device=device, region=region, name="NetDoping")) + is_depleted_0 = np.abs(electrons_0 - holes_0) < 0.1 * np.abs(doping_0) + dep_edges_0 = x_devsim[is_depleted_0] + if len(dep_edges_0) > 0 and np.min(dep_edges_0) <= 0.05: + v_punchthrough = 0.0 + else: + v_punchthrough = None # Calculate current density at 0V jn_edge_0 = np.array(devsim.get_edge_model_values(device=device, region=region, name=hf_opts['Jn'])) @@ -362,8 +441,8 @@ def build_and_solve_1d( else: step_size = abs_v / 20.0 - # For avalanche sweep, if current starts to rise, limit step size to 2V - if enable_avalanche and len(j_history) > 0 and j_history[-1] > 1e-4: + # For avalanche or BTBT sweep, if current starts to rise, limit step size to 2V + if (enable_avalanche or enable_btbt) and len(j_history) > 0 and j_history[-1] > 1e-4: step_size = min(step_size, 2.0) # Apply adaptive multiplier (from convergence failure backtracking) @@ -394,6 +473,16 @@ def build_and_solve_1d( v_history.append(v_current) j_history.append(j_val) + # Check for punch-through in sweep + if v_punchthrough is None: + electrons_sweep = np.array(devsim.get_node_model_values(device=device, region=region, name="Electrons")) + holes_sweep = np.array(devsim.get_node_model_values(device=device, region=region, name="Holes")) + doping_sweep = np.array(devsim.get_node_model_values(device=device, region=region, name="NetDoping")) + is_depleted_sweep = np.abs(electrons_sweep - holes_sweep) < 0.1 * np.abs(doping_sweep) + dep_edges_sweep = x_devsim[is_depleted_sweep] + if len(dep_edges_sweep) > 0 and np.min(dep_edges_sweep) <= 0.05: + v_punchthrough = v_current + # Check 1mA current limit current_ma = j_val * area_cm2 * 1000.0 if current_ma > 1.0: @@ -430,13 +519,22 @@ def build_and_solve_1d( jtot_edge = jn_edge + jp_edge g_av_edge = np.zeros_like(efield_edge) - if enable_avalanche: - try: - g_av_edge = np.array(devsim.get_edge_model_values(device=device, region=region, name="AvalancheGeneration")) - q = 1.6e-19 - g_av_edge = np.abs(g_av_edge) / q - except Exception: - pass + if enable_avalanche or enable_btbt: + q = 1.6e-19 + g_total = np.zeros_like(efield_edge) + if enable_avalanche: + try: + g_av = np.array(devsim.get_edge_model_values(device=device, region=region, name="AvalancheGeneration")) + g_total += np.abs(g_av) / q + except Exception: + pass + if enable_btbt: + try: + g_btbt = np.array(devsim.get_edge_model_values(device=device, region=region, name="BTBTGeneration")) + g_total += np.abs(g_btbt) / q + except Exception: + pass + g_av_edge = g_total # Final current density from right to left j_avg_rl = -np.mean(jtot_edge) @@ -476,7 +574,8 @@ def build_and_solve_1d( # Sweep history "v_history": v_history, "j_history": j_history, - "step_profiles": step_profiles + "step_profiles": step_profiles, + "v_punchthrough": v_punchthrough } @@ -501,6 +600,7 @@ if __name__ == "__main__": length=args['length'], process_steps=args['process_steps'], enable_avalanche=args['enable_avalanche'], + enable_btbt=args.get('enable_btbt', False), area_cm2=args['area_cm2'] ) diff --git a/gui1d/test_run.py b/gui1d/test_run.py index 87fec2c..1ae5988 100644 --- a/gui1d/test_run.py +++ b/gui1d/test_run.py @@ -1,6 +1,8 @@ # gui1d/test_run.py from solve_1d import build_and_solve_1d -print("Running test solve at V = -50.0V (Avalanche ON) with process steps...") +import numpy as np + +print("Running test sweep comparison up to 300V...") steps = [ { @@ -17,36 +19,46 @@ steps = [ ] try: - print("Call 1...") - res = build_and_solve_1d( - bias_target=50.0, + print("\n--- Solving case 1: Basic (Drift-Diffusion only) ---") + res_basic = build_and_solve_1d( + bias_target=300.0, substrate_type='n', substrate_doping=1e14, length=30.0, process_steps=steps, - enable_avalanche=True + enable_avalanche=False, + enable_btbt=False ) - print("Call 2...") - res2 = build_and_solve_1d( - bias_target=50.0, + print(f"Basic Solved Bias: {res_basic['v_solved']} V, Current Density: {res_basic['current_density']:.4e} A/cm^2") + + print("\n--- Solving case 2: Avalanche-only ---") + res_av = build_and_solve_1d( + bias_target=300.0, substrate_type='n', substrate_doping=1e14, length=30.0, process_steps=steps, - enable_avalanche=False + enable_avalanche=True, + enable_btbt=False ) - print("Call 3...") - res3 = build_and_solve_1d( - bias_target=5.0, + print(f"Avalanche Solved Bias: {res_av['v_solved']} V, Current Density: {res_av['current_density']:.4e} A/cm^2") + + print("\n--- Solving case 3: BTBT-only ---") + res_btbt = build_and_solve_1d( + bias_target=300.0, substrate_type='n', substrate_doping=1e14, length=30.0, process_steps=steps, - enable_avalanche=False + enable_avalanche=False, + enable_btbt=True ) - print("All 3 calls passed successfully!") + print(f"BTBT Solved Bias: {res_btbt['v_solved']} V, Current Density: {res_btbt['current_density']:.4e} A/cm^2") + + print("\nAll comparison runs finished successfully!") except Exception as e: print("Test failed with error:") import traceback traceback.print_exc() + diff --git a/make_scripts.py b/make_scripts.py index 763a380..5b1c6d9 100644 --- a/make_scripts.py +++ b/make_scripts.py @@ -28,7 +28,7 @@ recon_logic = """ if v_current >= next_recon_v: state_data = save_state(device) seed_data = {"voltage": v_current, "step_size": step_size, "state": state_data} - seed_filename = f"seed_{int(next_recon_v)}V.pkl" + seed_filename = f"{OUT_DIR}seed_{int(next_recon_v)}V.pkl" with open(seed_filename, "wb") as f: pickle.dump(seed_data, f) print(f"\\n--- RECON PROBE at {v_current:.2f} V ---") @@ -58,15 +58,15 @@ recon_logic = """ ia_p_p12 = devsim.get_contact_current(device=device, contact="MT1_P12_Si", equation="HoleContinuityEquation") av_curr = ia_n_si + ia_p_si + ia_n_p12 + ia_p_p12 print(f"Avalanche Current at {v_current:.2f} V: {av_curr:.4e} A") - with open("recon_avalanche.log", "a") as f: + with open(f"{OUT_DIR}recon_avalanche.log", "a") as f: f.write(f"{v_current:.2f}\\t{av_curr:.4e}\\n") else: print("Avalanche failed to converge.") - with open("recon_avalanche.log", "a") as f: + with open(f"{OUT_DIR}recon_avalanche.log", "a") as f: f.write(f"{v_current:.2f}\\tFAILED\\n") except devsim.error: print("Avalanche failed to converge.") - with open("recon_avalanche.log", "a") as f: + with open(f"{OUT_DIR}recon_avalanche.log", "a") as f: f.write(f"{v_current:.2f}\\tFAILED\\n") # Restore state and Turn OFF Avalanche @@ -90,7 +90,7 @@ current_list = [0.0] # Recon variables next_recon_v = 50.0 -with open("recon_avalanche.log", "w") as f: +with open(f"{OUT_DIR}recon_avalanche.log", "w") as f: f.write("Voltage(V)\\tAvalancheCurrent(A)\\n") """ diff --git a/physics/new_physics.py b/physics/new_physics.py index b25132e..fa95472 100644 --- a/physics/new_physics.py +++ b/physics/new_physics.py @@ -527,3 +527,34 @@ def CreateAvalancheGeneration(device, region, Jn, Jp): CreateEdgeModelDerivatives(device, region, "AvalancheGeneration_p", eq_p, i) return "AvalancheGeneration" + + +def CreateBTBTGeneration(device, region): + ''' + Band-to-Band Tunneling (BTBT) model using Hurkx formulation. + Generates edge models `BTBTGeneration` and `BTBTGeneration_p` to be used + in continuity equations. + ''' + # Set Hurkx parameters + set_parameter(device=device, region=region, name="A_btbt", value=4.0e14) # cm^-1 s^-1 V^-2 + set_parameter(device=device, region=region, name="B_btbt", value=1.9e7) # V/cm + + # Calculate Electric Field Magnitude if not already created + if not InEdgeModelList(device, region, "E_mag"): + CreateEdgeModel(device, region, "E_mag", "(EField^2 + 1e-4)^0.5") + for i in ("Potential",): + CreateEdgeModelDerivatives(device, region, "E_mag", "(EField^2 + 1e-4)^0.5", i) + + # Calculate BTBT Generation rate (multiplied by q, so unit is A/cm^3) + eq_n = "q * A_btbt * (E_mag^2) * exp(-B_btbt / E_mag)" + CreateEdgeModel(device, region, "BTBTGeneration", eq_n) + + eq_p = "-q * A_btbt * (E_mag^2) * exp(-B_btbt / E_mag)" + CreateEdgeModel(device, region, "BTBTGeneration_p", eq_p) + + for i in ("Potential",): + CreateEdgeModelDerivatives(device, region, "BTBTGeneration", eq_n, i) + CreateEdgeModelDerivatives(device, region, "BTBTGeneration_p", eq_p, i) + + return "BTBTGeneration" + diff --git a/preview_doping_2d.py b/preview_doping_2d.py index ae75716..44f68d0 100644 --- a/preview_doping_2d.py +++ b/preview_doping_2d.py @@ -2,12 +2,17 @@ import devsim import numpy as np import matplotlib.pyplot as plt import matplotlib.tri as tri +import os +import sys +DEV_DIR = os.environ.get("DEV_DIR", "devices/Triac_rp") +sys.path.insert(0, os.path.abspath(DEV_DIR)) from device_config import * device = "device_2d" # 1. Load the mesh -devsim.create_gmsh_mesh(mesh=device, file="device_2d.msh") +mesh_file = os.path.join(DEV_DIR, "device_2d.msh") +devsim.create_gmsh_mesh(mesh=device, file=mesh_file) devsim.add_gmsh_region(mesh=device, gmsh_name="Silicon", region="Silicon", material="Silicon") devsim.add_gmsh_region(mesh=device, gmsh_name="Oxide", region="Oxide", material="Oxide") devsim.add_gmsh_region(mesh=device, gmsh_name="Molding", region="Molding", material="Molding") @@ -41,58 +46,63 @@ devsim.finalize_mesh(mesh=device) devsim.create_device(mesh=device, device=device) # 2. Define Doping Profiles using sub-models to avoid long strings -# Substrate (N-type) -devsim.node_model(device=device, region="Silicon", name="nD_sub", equation=f"{N_SUB}") - -# Helper to generate 2D erfc profile string -def get_erfc_expr(peak, x1, x2, hdiff, vdiff): - return f"{peak} * erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))" - -# P-well profiles (p11, p12, p13 on both sides) -# p11 -p11_left_expr = get_erfc_expr(P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p11_right_expr = get_erfc_expr(P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p11_l", equation=p11_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p11_r", equation=p11_right_expr) - -# p12 -p12_left_expr = get_erfc_expr(P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p12_right_expr = get_erfc_expr(P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p12_l", equation=p12_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p12_r", equation=p12_right_expr) - -# p13 -p13_left_expr = get_erfc_expr(P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p13_right_expr = get_erfc_expr(P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p13_l", equation=p13_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p13_r", equation=p13_right_expr) - -# N+ profiles -nplus_left_expr = get_erfc_expr(NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) -nplus_right_expr = get_erfc_expr(NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nD_nplus_l", equation=nplus_left_expr) -devsim.node_model(device=device, region="Silicon", name="nD_nplus_r", equation=nplus_right_expr) - -# MRING N+ profiles -mring_l_expr = get_erfc_expr(NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) -mring_r_expr = get_erfc_expr(NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nD_mring_l", equation=mring_l_expr) -devsim.node_model(device=device, region="Silicon", name="nD_mring_r", equation=mring_r_expr) - -# Combine into Donors and Acceptors -devsim.node_model(device=device, region="Silicon", name="Donors", - equation="nD_sub + nD_nplus_l + nD_nplus_r + nD_mring_l + nD_mring_r") -devsim.node_model(device=device, region="Silicon", name="Acceptors", - equation="1e10 + nA_p11_l + nA_p11_r + nA_p12_l + nA_p12_r + nA_p13_l + nA_p13_r") - -# NetDoping -devsim.node_model(device=device, region="Silicon", name="NetDoping", equation="Donors - Acceptors") -devsim.node_model(device=device, region="Silicon", name="LogNetDoping", equation="asinh(NetDoping / 2.0) / log(10.0)") -devsim.node_model(device=device, region="Silicon", name="LogAcceptors", equation="log(Acceptors) / log(10.0)") +if os.environ.get("USE_PCAD", "false").lower() == "true": + from device_pcad_config import apply_pcad_doping_2d + apply_pcad_doping_2d(device, region="Silicon") + devsim.node_model(device=device, region="Silicon", name="LogAcceptors", equation="log(Acceptors) / log(10.0)") +else: + # Substrate (N-type) + devsim.node_model(device=device, region="Silicon", name="nD_sub", equation=f"{N_SUB}") + + # Helper to generate 2D erfc profile string + def get_erfc_expr(peak, x1, x2, hdiff, vdiff): + return f"{peak} * erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))" + + # P-well profiles (p11, p12, p13 on both sides) + # p11 + p11_left_expr = get_erfc_expr(P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p11_right_expr = get_erfc_expr(P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p11_l", equation=p11_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p11_r", equation=p11_right_expr) + + # p12 + p12_left_expr = get_erfc_expr(P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p12_right_expr = get_erfc_expr(P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p12_l", equation=p12_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p12_r", equation=p12_right_expr) + + # p13 + p13_left_expr = get_erfc_expr(P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p13_right_expr = get_erfc_expr(P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p13_l", equation=p13_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p13_r", equation=p13_right_expr) + + # N+ profiles + nplus_left_expr = get_erfc_expr(NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) + nplus_right_expr = get_erfc_expr(NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nD_nplus_l", equation=nplus_left_expr) + devsim.node_model(device=device, region="Silicon", name="nD_nplus_r", equation=nplus_right_expr) + + # MRING N+ profiles + mring_l_expr = get_erfc_expr(NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) + mring_r_expr = get_erfc_expr(NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nD_mring_l", equation=mring_l_expr) + devsim.node_model(device=device, region="Silicon", name="nD_mring_r", equation=mring_r_expr) + + # Combine into Donors and Acceptors + devsim.node_model(device=device, region="Silicon", name="Donors", + equation="nD_sub + nD_nplus_l + nD_nplus_r + nD_mring_l + nD_mring_r") + devsim.node_model(device=device, region="Silicon", name="Acceptors", + equation="1e10 + nA_p11_l + nA_p11_r + nA_p12_l + nA_p12_r + nA_p13_l + nA_p13_r") + + # NetDoping + devsim.node_model(device=device, region="Silicon", name="NetDoping", equation="Donors - Acceptors") + devsim.node_model(device=device, region="Silicon", name="LogNetDoping", equation="asinh(NetDoping / 2.0) / log(10.0)") + devsim.node_model(device=device, region="Silicon", name="LogAcceptors", equation="log(Acceptors) / log(10.0)") # Write Tecplot output for ParaView -devsim.write_devices(file="device_2d.tec", type="tecplot") -devsim.write_devices(file="preview.tec", type="tecplot") +devsim.write_devices(file=os.path.join(DEV_DIR, "device_2d.tec"), type="tecplot") +devsim.write_devices(file=os.path.join(DEV_DIR, "preview.tec"), type="tecplot") print("Saved device_2d.tec and preview.tec") # 4. Generate a 2D Plot with Matplotlib to verify the doping profile @@ -202,6 +212,7 @@ ax2.set_xlim(-W_SIM / um, W_SIM / um) ax2.set_ylim(H_SI/um + 15.0, -110.0) plt.tight_layout() -plt.savefig('doping_2d.png', dpi=300) +png_out = os.path.join(DEV_DIR, "doping_2d.png") +plt.savefig(png_out, dpi=300) plt.close() -print("Plot saved to doping_2d.png") +print(f"Plot saved to {png_out}") diff --git a/project_discussion.md b/project_discussion.md index c137fdf..1c7ade5 100644 --- a/project_discussion.md +++ b/project_discussion.md @@ -1365,4 +1365,71 @@ $$\Delta V < W_{gb} \cdot \sqrt{\frac{2 q \cdot N_{SUB} \cdot V}{\epsilon}}$$ - 這徹底解決了高反向偏壓下雪崩電流倒灌變負、以及在大偏壓網格細化時的發散問題，使 1D 和 2D 雪崩模擬皆能 100% 收斂且符合物理實相。 +### 25. 關於二極體與 BJT 穿通機制之物理本質差異，以及未來 Zener 穿隧效應（BTBT）之擴充建議（2026-06-20） +在進行 1D Diode 模擬與高壓穿通分析時，我們深入討論了二極體與三極體（BJT）在穿通（Punch-through）狀態下的電學行為與數值模擬特徵，並為未來引入低壓齊納穿隧效應提出了具體的架構性建議。 + +#### 25.1 穿通（Punch-through）與擊穿（Breakdown）的物理本質差異 + +1. **二極體 ($P^+ - N - N^+$) 穿通的特徵**： + * **物理機制**：當反向偏壓增加時，空乏區延伸至同類載流子接觸區（例如基底 $N^-$ 空乏區觸及底部的 $N^+$ 歐姆區，或 P-well 被吹透至 $x=0$ 的陽極金屬）。此時，元件內部的電場分布會由「三角形」轉變為「梯形」。 + * **電流行為**：在無雪崩效應的理想 Drift-Diffusion 框架下，穿通本身**不會**導致漏電流激增。因為理想歐姆接觸的邊界條件強制鎖定了少數載流子濃度（例如 $n(0) \approx 10^3 \text{ cm}^{-3}$），金屬接觸面無法源源不絕地向空乏區注入電子。因此，電流依然只由微弱的熱產生控制，維持在極低且飽和的狀態。這與真實物理中，二極體穿通後需要繼續升高電壓觸發**雪崩擊穿（Avalanche Breakdown）**電流才會暴增的現象完全吻合。 + +2. **BJT ($N^+ - P - N$) 穿通的特徵**： + * **物理機制**：在三層結構中，穿通空間意味着中間阻擋層（Base, P 區）被完全空乏，使得發射極（Emitter, $N^+$）與集極（Collector, $N$）之間的勢壘（Barrier）徹底坍塌。 + * **電流行為**：由於發射極接觸面為 $N^+$，邊界上存在極高濃度的電子源（$n(0) \approx 10^{19} \text{ cm}^{-3}$）。一旦勢壘消除，大量電子會瞬間注入並被電場掃向 Collector，導致 I-V 曲線在穿通電壓處發生**指數級的電流暴增**（即使沒有開啟雪崩效應）。 + +#### 25.2 目前已實作之物理機制彙整 + +當前專案中已實作且驗證的物理機制包括： +* **Scharfetter-Gummel 漂移-擴散傳導**（高電場數值穩定） +* **Arora 低場與 Canali 高場飽和遷移率模型**（速度飽和） +* **Shockley-Read-Hall (SRH) 熱複合-產生模型** +* **Chynoweth 碰撞游離雪崩模型**（電荷守恆對稱宣告） + +#### 25.3 未來可擴充之機制：Zener 穿隧效應 (BTBT) 建議 + +當 TVS 元件設計於低壓工作區（如 $V_Z < 5\text{ V}$）時，接面兩側摻雜極高，空乏區極窄，會觸發量子力學的**能帶間穿隧（Band-to-Band Tunneling, BTBT）**。目前專案尚未包含此機制，未來可依循以下方案於 `physics/new_physics.py` 中掛載： + +1. **宣告 Hurkx 穿隧模型 (Hurkx BTBT Model)**： + 利用局部電場強度 $E$，定義 BTBT 產生率模型： + $$G_{BTBT} = A \cdot E^2 \cdot \exp\left( - \frac{B}{E} \right)$$ + 在 Python 中宣告對應的邊緣模型（Edge Model）與對 `Potential` 的 Jacobian 偏微分導數。 + +2. **將穿隧源項納入連續方程式**： + 在 `ElectronContinuityEquation` 與 `HoleContinuityEquation` 中，將 `BTBT` 產生率累加至 `edge_volume_model` 之中（與原有的 `AvalancheGeneration` 疊加）： + * 電子方程：`edge_volume_model = "AvalancheGeneration + BTBTGeneration"` + * 電洞方程：`edge_volume_model = "AvalancheGeneration_p - BTBTGeneration"` (維持電荷守恆) + +這將使元件在低壓重摻雜下能正確模擬出穿隧漏電與穩壓特性，是專案後續向低壓 TVS 元件擴展的關鍵物理拼圖。 + + +### 26. 關於未來 2D 製程配置與元件建模之設計討論（2026-06-25） + +為了在將來導入 2D 模擬時，能夠更靈活、更貼近實際製程思維地建立元件摻雜配置（Doping Configuration），我們討論並整理了以下 2D 元件設定的需求作為後續的開發計畫 (To-Do List)。 + +#### 26.1 2D 元件配置之設計思維與需求 +未來在設計 2D `device_config.py` 與對應的圖形介面時，製程建模的前置作業將以 **「2D 幾何與局部開窗」** 的思維來獨立設計，暫不強求與 1D 模擬邏輯進行直接的連動。 + +具體製程步驟之配置應支援以下三種核心類型： + +1. **局部摻雜步驟 (Selective Doping Step: Implant / Diffusion)**： + 每一步驟皆可包含以下參數： + * **摻雜類型 (Type)**：`n` 或 `p`。 + * **方法 (Method)**：`Implant` (離子佈植) 或 `Diffusion` (熱擴散)。 + * **劑量/濃度 (Dosage / Surface Cs)**：定量摻雜參數。 + * **能量 (Energy)**：離子佈植能量（若為 Implant）。 + * **熱預算 (Thermal Budget)**：溫度（Temperature）與時間（Time）。 + * **局部開窗區域 (Opening Areas)**：指定水平方向的多個開窗區間，例如 `x1 ~ x2`、`x3 ~ x4` 等，只有在開窗區內才會引入雜質。 + +2. **局部埋層步驟 (Buried Layer Step)**： + 用以直接在 Epi 下方或元件特定深處定義局部高濃度埋層。參數應包含： + * **埋層物理參數**：峰值濃度、縱向深度、縱向與橫向擴散長度等。 + * **影響區域 (Affected Areas)**：水平方向的作用區間，例如 `x11 ~ x12`、`x13 ~ x14` 等。 + +3. **磊晶成長步驟 (Epi Growth Step)**： + * **磊晶規格**：厚度、摻雜類型與濃度等。 + * 此步驟通常是全域性的，用於覆蓋先前已進行局部摻雜/埋層的表面，並在後續熱處理中做為雜質雙向擴散的介質。 + +#### 26.2 實現路徑之評估 +在 2D 實現此設計時，我們將繼續秉持**不直接在 DEVSIM 中執行網格形變/動態擴散偏微分求解**的原則，以降低實作複雜度。取而代之的是，在靜態的 2D 網格中，利用解析函數（Analytical Functions）結合水平方向的 `erf` 邊緣過渡函數，動態拼接出符合各步驟開窗範圍的 2D 摻雜分佈。 diff --git a/resume_run.py b/resume_run.py index eb981d6..5649ce7 100644 --- a/resume_run.py +++ b/resume_run.py @@ -4,6 +4,10 @@ import sys import glob import gc +DEV_DIR = os.environ.get("DEV_DIR", "devices/Triac_rp") +sys.path.insert(0, os.path.abspath(DEV_DIR)) +OUT_DIR = os.path.join(os.environ.get("OUT_DIR", os.path.join(DEV_DIR, "output_this_run")), "") + # Limit the thread count for parallel solvers to prevent WSL from resource starvation/disconnecting os.environ["OMP_NUM_THREADS"] = "6" os.environ["MKL_NUM_THREADS"] = "6" @@ -22,12 +26,13 @@ import devsim import numpy as np is_avalanche_enabled = os.environ.get("AVALANCHE", "false").lower() == "true" +is_btbt_enabled = os.environ.get("BTBT", "false").lower() == "true" refine_v_step = float(os.environ.get("REFINE_V_STEP", "50.0")) is_refine_enabled = os.environ.get("REFINE", "false").lower() == "true" if refine_v_step < 1.0: is_refine_enabled = False print(f"refine_v_step < 1.0 V detected ({refine_v_step} V): Dynamic refinement is completely DISABLED.") -print(f"Option: AVALANCHE={is_avalanche_enabled}, REFINE={is_refine_enabled}, REFINE_V_STEP={refine_v_step}V") +print(f"Option: AVALANCHE={is_avalanche_enabled}, BTBT={is_btbt_enabled}, REFINE={is_refine_enabled}, REFINE_V_STEP={refine_v_step}V") import matplotlib.pyplot as plt import time @@ -59,7 +64,7 @@ if specified_checkpoint: print(f"Error: Specified checkpoint path '{specified_checkpoint}' does not exist.") else: # 2. 自動搜尋邏輯（原本的機制） - search_dirs = ["output_this_run/"] + search_dirs = [OUT_DIR] for d in glob.glob("output_*"): if os.path.isdir(d): search_dirs.append(d + "/") @@ -102,18 +107,18 @@ latest_mtime, latest_voltage, latest_cp, cp_data = valid_checkpoints[0] cp_dir = os.path.dirname(latest_cp) if cp_dir: - OUT_DIR = cp_dir + "/" -else: - OUT_DIR = "output_this_run/" + OUT_DIR = os.path.join(cp_dir, "") print(f"Resuming outputs to directory: {OUT_DIR}") os.makedirs(OUT_DIR, exist_ok=True) import shutil try: - if os.path.exists("device_2d.msh"): - shutil.copy2("device_2d.msh", os.path.join(OUT_DIR, "device_2d.msh")) - if os.path.exists("device_config.py"): - shutil.copy2("device_config.py", os.path.join(OUT_DIR, "device_config.py")) + mesh_orig = os.path.join(DEV_DIR, "device_2d.msh") + if os.path.exists(mesh_orig): + shutil.copy2(mesh_orig, os.path.join(OUT_DIR, "device_2d.msh")) + config_orig = os.path.join(DEV_DIR, "device_config.py") + if os.path.exists(config_orig): + shutil.copy2(config_orig, os.path.join(OUT_DIR, "device_config.py")) if __file__: shutil.copy2(__file__, os.path.join(OUT_DIR, os.path.basename(__file__))) except Exception as e: @@ -135,7 +140,7 @@ silicon_state_len = len(cp_data["state"]["Silicon"]["Potential"]) try: # Load coarse mesh temporarily to count Silicon nodes - devsim.create_gmsh_mesh(mesh="temp_check", file="device_2d.msh") + devsim.create_gmsh_mesh(mesh="temp_check", file=os.path.join(DEV_DIR, "device_2d.msh")) devsim.add_gmsh_region(mesh="temp_check", gmsh_name="Silicon", region="Silicon", material="Silicon") devsim.add_gmsh_region(mesh="temp_check", gmsh_name="Oxide", region="Oxide", material="Oxide") devsim.add_gmsh_region(mesh="temp_check", gmsh_name="Molding", region="Molding", material="Molding") @@ -160,7 +165,7 @@ if is_refined: cp_dir = os.path.dirname(latest_cp) mesh_candidates = [ os.path.join(cp_dir, "device_2d_refined.msh") if cp_dir else None, - os.path.join("output_this_run", "device_2d_refined.msh") + os.path.join(OUT_DIR, "device_2d_refined.msh") ] mesh_file = None for cand in mesh_candidates: @@ -176,10 +181,10 @@ if is_refined: raise RuntimeError("Refined checkpoint requires device_2d_refined.msh, but it was not found.") print(f"Refined checkpoint detected. Loading refined mesh: {mesh_file}") else: - mesh_file = "device_2d.msh" + mesh_file = os.path.join(DEV_DIR, "device_2d.msh") print(f"Standard checkpoint detected. Loading base mesh: {mesh_file}") -sys.path.append("/home/pchan/devsim2026") +# sys.path already configured at top from device_config import * from physics.model_create import * from physics.new_physics import * @@ -216,42 +221,46 @@ devsim.finalize_mesh(mesh=device) devsim.create_device(mesh=device, device=device) # 2. Set up doping in Silicon region -devsim.node_model(device=device, region="Silicon", name="nD_sub", equation=f"{N_SUB}") +if os.environ.get("USE_PCAD", "false").lower() == "true": + from device_pcad_config import apply_pcad_doping_2d + apply_pcad_doping_2d(device, region="Silicon") +else: + devsim.node_model(device=device, region="Silicon", name="nD_sub", equation=f"{N_SUB}") -def get_erfc_expr(peak, x1, x2, hdiff, vdiff): - return f"{peak} * erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))" + def get_erfc_expr(peak, x1, x2, hdiff, vdiff): + return f"{peak} * erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))" -p11_left_expr = get_erfc_expr(P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p11_right_expr = get_erfc_expr(P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p11_l", equation=p11_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p11_r", equation=p11_right_expr) + p11_left_expr = get_erfc_expr(P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p11_right_expr = get_erfc_expr(P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p11_l", equation=p11_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p11_r", equation=p11_right_expr) -p12_left_expr = get_erfc_expr(P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p12_right_expr = get_erfc_expr(P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p12_l", equation=p12_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p12_r", equation=p12_right_expr) + p12_left_expr = get_erfc_expr(P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p12_right_expr = get_erfc_expr(P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p12_l", equation=p12_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p12_r", equation=p12_right_expr) -p13_left_expr = get_erfc_expr(P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p13_right_expr = get_erfc_expr(P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p13_l", equation=p13_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p13_r", equation=p13_right_expr) + p13_left_expr = get_erfc_expr(P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p13_right_expr = get_erfc_expr(P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p13_l", equation=p13_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p13_r", equation=p13_right_expr) -nplus_left_expr = get_erfc_expr(NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) -nplus_right_expr = get_erfc_expr(NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nD_nplus_l", equation=nplus_left_expr) -devsim.node_model(device=device, region="Silicon", name="nD_nplus_r", equation=nplus_right_expr) + nplus_left_expr = get_erfc_expr(NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) + nplus_right_expr = get_erfc_expr(NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nD_nplus_l", equation=nplus_left_expr) + devsim.node_model(device=device, region="Silicon", name="nD_nplus_r", equation=nplus_right_expr) -mring_l_expr = get_erfc_expr(NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) -mring_r_expr = get_erfc_expr(NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nD_mring_l", equation=mring_l_expr) -devsim.node_model(device=device, region="Silicon", name="nD_mring_r", equation=mring_r_expr) + mring_l_expr = get_erfc_expr(NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) + mring_r_expr = get_erfc_expr(NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nD_mring_l", equation=mring_l_expr) + devsim.node_model(device=device, region="Silicon", name="nD_mring_r", equation=mring_r_expr) -devsim.node_model(device=device, region="Silicon", name="Donors", - equation="nD_sub + nD_nplus_l + nD_nplus_r + nD_mring_l + nD_mring_r") -devsim.node_model(device=device, region="Silicon", name="Acceptors", - equation="1e10 + nA_p11_l + nA_p11_r + nA_p12_l + nA_p12_r + nA_p13_l + nA_p13_r") -devsim.node_model(device=device, region="Silicon", name="NetDoping", equation="Donors - Acceptors") -devsim.node_model(device=device, region="Silicon", name="LogNetDoping", equation="asinh(NetDoping / 2.0) / log(10.0)") + devsim.node_model(device=device, region="Silicon", name="Donors", + equation="nD_sub + nD_nplus_l + nD_nplus_r + nD_mring_l + nD_mring_r") + devsim.node_model(device=device, region="Silicon", name="Acceptors", + equation="1e10 + nA_p11_l + nA_p11_r + nA_p12_l + nA_p12_r + nA_p13_l + nA_p13_r") + devsim.node_model(device=device, region="Silicon", name="NetDoping", equation="Donors - Acceptors") + devsim.node_model(device=device, region="Silicon", name="LogNetDoping", equation="asinh(NetDoping / 2.0) / log(10.0)") # 3. Initialize electrostatic potential simulation (Poisson only) CreateSolution(device, "Silicon", "Potential") @@ -384,17 +393,33 @@ print("Skipping initial Drift-Diffusion solve, preparing to restore from checkpo # Switch continuity and potential equations for the bias sweep print("Configuring continuity and potential equations for the bias sweep (min_error=1e5, positive update)...") -# Instantiate Avalanche (Impact Ionization) edge generation model if enabled -if is_avalanche_enabled: - CreateAvalancheGeneration(device, "Silicon", opts['Jn'], opts['Jp']) +# Instantiate Avalanche and BTBT generation models +CreateAvalancheGeneration(device, "Silicon", opts['Jn'], opts['Jp']) +CreateBTBTGeneration(device, "Silicon") av_model_n = "AvalancheGeneration" if is_avalanche_enabled else "" av_model_p = "AvalancheGeneration_p" if is_avalanche_enabled else "" +btbt_model_n = "BTBTGeneration" if is_btbt_enabled else "" +btbt_model_p = "BTBTGeneration_p" if is_btbt_enabled else "" + +if av_model_n and btbt_model_n: + main_edge_volume_model_n = f"{av_model_n} + {btbt_model_n}" + main_edge_volume_model_p = f"{av_model_p} + {btbt_model_p}" +elif av_model_n: + main_edge_volume_model_n = av_model_n + main_edge_volume_model_p = av_model_p +elif btbt_model_n: + main_edge_volume_model_n = btbt_model_n + main_edge_volume_model_p = btbt_model_p +else: + main_edge_volume_model_n = "" + main_edge_volume_model_p = "" + devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons", - time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=av_model_n, + time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=main_edge_volume_model_n, variable_update="positive", node_model="ElectronGeneration", min_error=1e5) devsim.equation(device=device, region="Silicon", name="HoleContinuityEquation", variable_name="Holes", - time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=av_model_p, + time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=main_edge_volume_model_p, variable_update="positive", node_model="HoleGeneration", min_error=1e5) devsim.equation(device=device, region="Silicon", name="PotentialEquation", variable_name="Potential", node_model="PotentialNodeCharge", edge_model="DField", variable_update="default", min_error=1e-3) @@ -445,7 +470,7 @@ print(f"Resuming from V = {v_current:.4f} V, step = {step_size:.4f} V") # 為了讓 log 連續，若 checkpoint 來自其他目錄，將舊日誌複製至當前的 output_this_run/ 中 cp_dir = os.path.dirname(latest_cp) if cp_dir and cp_dir != OUT_DIR.rstrip("/"): - for log_name in ["simulation_time.log", "recon_avalanche.log"]: + for log_name in ["simulation_time.log", "recon_avalanche.log", "recon_btbt.log", "recon_av_btbt.log"]: old_log = os.path.join(cp_dir, log_name) new_log = os.path.join(OUT_DIR, log_name) if os.path.exists(old_log) and not os.path.exists(new_log): @@ -464,9 +489,9 @@ for c in ["MT1_Si", "MT1_P12_Si", "MT1_Ox", "MT1_Mold"]: restore_state(device, state) # File logging setup (append mode for resume) -time_log = open(f"{OUT_DIR}simulation_time.log", "a", buffering=1) -with open(f"{OUT_DIR}recon_avalanche.log", "a") as f: - pass # Append mode for recon +for log_name in ["recon_avalanche.log", "recon_btbt.log", "recon_av_btbt.log"]: + with open(f"{OUT_DIR}{log_name}", "a") as f: + pass start_sweep_time = time.time() @@ -496,6 +521,9 @@ saved_current_targets = {t for t in TEC_CURRENT_TARGETS if any(abs(i) >= t for i seed_save_targets = [5.0, 25.0, 45.0, 95.0, 195.0, 395.0, 595.0, 795.0, 995.0, 1195.0] saved_seeds = {t for t in seed_save_targets if v_current >= t} +# Open time log +time_log = open(f"{OUT_DIR}simulation_time.log", "a", buffering=1) + while v_current < v_target: v_next = min(v_current + step_size, v_target) @@ -510,13 +538,11 @@ while v_current < v_target: iters1 = 15 try: # Always use log_damp preconditioning for Electron/Hole continuity equations in Stage 1 - av_model_n = "AvalancheGeneration" if is_avalanche_enabled else "" - av_model_p = "AvalancheGeneration_p" if is_avalanche_enabled else "" devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons", - time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=av_model_n, + time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=main_edge_volume_model_n, variable_update="log_damp", node_model="ElectronGeneration", min_error=1e5) devsim.equation(device=device, region="Silicon", name="HoleContinuityEquation", variable_name="Holes", - time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=av_model_p, + time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=main_edge_volume_model_p, variable_update="log_damp", node_model="HoleGeneration", min_error=1e5) res1 = devsim.solve(type="dc", absolute_error=1e10, relative_error=1e-3, charge_error=1e12, maximum_iterations=5, rollback=False, info=True) iters1 = len(res1.get("iterations", [])) @@ -524,13 +550,11 @@ while v_current < v_target: pass # Ignore non-convergence in pre-conditioning finally: # Always revert to positive variable update for Stage 2 precision Newton - av_model_n = "AvalancheGeneration" if is_avalanche_enabled else "" - av_model_p = "AvalancheGeneration_p" if is_avalanche_enabled else "" devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons", - time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=av_model_n, + time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=main_edge_volume_model_n, variable_update="positive", node_model="ElectronGeneration", min_error=1e5) devsim.equation(device=device, region="Silicon", name="HoleContinuityEquation", variable_name="Holes", - time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=av_model_p, + time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=main_edge_volume_model_p, variable_update="positive", node_model="HoleGeneration", min_error=1e5) import psutil @@ -625,7 +649,7 @@ while v_current < v_target: import dynamic_refine try: # 1. 執行網格自適應重劃與狀態插值 - refined_device, refined_opts = dynamic_refine.refine_and_interpolate(device, v_current, is_avalanche_enabled=is_avalanche_enabled, time_log=time_log, out_dir=OUT_DIR) + refined_device, refined_opts = dynamic_refine.refine_and_interpolate(device, v_current, is_avalanche_enabled=is_avalanche_enabled, is_btbt_enabled=is_btbt_enabled, time_log=time_log, out_dir=OUT_DIR) device = refined_device opts = refined_opts just_refined = True @@ -640,15 +664,19 @@ while v_current < v_target: print(f"\n--- RECON PROBE at {v_current:.2f} V ---") print(f"Saved refined seed to {seed_filename}") - if is_avalanche_enabled: - # Turn ON Avalanche (Use log_damp for Stage 1 pre-conditioning stability) + # Helper to run a recon probe + def run_recon_probe(probe_name, model_n, model_p, log_file): + print(f" Running Recon Probe ({probe_name})...") + # Temporarily set equations devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons", - time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model="AvalancheGeneration", - variable_update="log_damp", node_model="ElectronGeneration", min_error=1e5) + time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=model_n, + variable_update="positive", node_model="ElectronGeneration", min_error=1e5) devsim.equation(device=device, region="Silicon", name="HoleContinuityEquation", variable_name="Holes", - time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model="AvalancheGeneration_p", - variable_update="log_damp", node_model="HoleGeneration", min_error=1e5) + time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=model_p, + variable_update="positive", node_model="HoleGeneration", min_error=1e5) + converged = False + total_curr = 0.0 try: # Stage 1 pre-conditioning try: @@ -662,10 +690,10 @@ while v_current < v_target: # Switch to positive update for Stage 2 precision Newton solve devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons", - time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model="AvalancheGeneration", + time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=model_n, variable_update="positive", node_model="ElectronGeneration", min_error=1e5) devsim.equation(device=device, region="Silicon", name="HoleContinuityEquation", variable_name="Holes", - time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model="AvalancheGeneration_p", + time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=model_p, variable_update="positive", node_model="HoleGeneration", min_error=1e5) # Stage 2 solve (Strict Tolerance) @@ -674,34 +702,40 @@ while v_current < v_target: mem_recon_stage2 = psutil.Process(os.getpid()).memory_info().rss / (1024**3) print(f"Recon Stage 2 Memory Usage: {mem_recon_stage2:.1f} GB") if res_av.get("converged", False): - # Measure Avalanche current ia_n_si = devsim.get_contact_current(device=device, contact="MT1_Si", equation="ElectronContinuityEquation") ia_p_si = devsim.get_contact_current(device=device, contact="MT1_Si", equation="HoleContinuityEquation") ia_n_p12 = devsim.get_contact_current(device=device, contact="MT1_P12_Si", equation="ElectronContinuityEquation") ia_p_p12 = devsim.get_contact_current(device=device, contact="MT1_P12_Si", equation="HoleContinuityEquation") - av_curr = ia_n_si + ia_p_si + ia_n_p12 + ia_p_p12 - print(f"Avalanche Current at {v_current:.2f} V: {av_curr:.4e} A") - with open(f"{OUT_DIR}recon_avalanche.log", "a") as f: - f.write(f"{v_current:.2f}\t{av_curr:.4e}\n") - else: - print("Avalanche failed to converge.") - with open(f"{OUT_DIR}recon_avalanche.log", "a") as f: - f.write(f"{v_current:.2f}\tFAILED\n") + total_curr = ia_n_si + ia_p_si + ia_n_p12 + ia_p_p12 + converged = True except devsim.error: - print("Avalanche failed to converge.") - with open(f"{OUT_DIR}recon_avalanche.log", "a") as f: - f.write(f"{v_current:.2f}\tFAILED\n") + pass - # Restore state and Reset Avalanche state to main loop configuration + # Restore state immediately after solve attempts restore_state(device, state) - av_model_n = "AvalancheGeneration" if is_avalanche_enabled else "" - av_model_p = "AvalancheGeneration_p" if is_avalanche_enabled else "" - devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons", - time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=av_model_n, - variable_update="positive", node_model="ElectronGeneration", min_error=1e5) - devsim.equation(device=device, region="Silicon", name="HoleContinuityEquation", variable_name="Holes", - time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=av_model_p, - variable_update="positive", node_model="HoleGeneration", min_error=1e5) + + with open(log_file, "a") as f: + if converged: + print(f" {probe_name} Current at {v_current:.2f} V: {total_curr:.4e} A") + f.write(f"{v_current:.2f}\t{total_curr:.4e}\n") + else: + print(f" {probe_name} failed to converge.") + f.write(f"{v_current:.2f}\tFAILED\n") + + # 1. Avalanche Only + run_recon_probe("Avalanche", "AvalancheGeneration", "AvalancheGeneration_p", f"{OUT_DIR}recon_avalanche.log") + # 2. BTBT Only + run_recon_probe("BTBT", "BTBTGeneration", "BTBTGeneration_p", f"{OUT_DIR}recon_btbt.log") + # 3. Avalanche + BTBT + run_recon_probe("Av+BTBT", "AvalancheGeneration + BTBTGeneration", "AvalancheGeneration_p + BTBTGeneration_p", f"{OUT_DIR}recon_av_btbt.log") + + # Revert equations to main sweep configuration + devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons", + time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model=main_edge_volume_model_n, + variable_update="positive", node_model="ElectronGeneration", min_error=1e5) + devsim.equation(device=device, region="Silicon", name="HoleContinuityEquation", variable_name="Holes", + time_node_model="PCharge", edge_model=opts['Jp'], edge_volume_model=main_edge_volume_model_p, + variable_update="positive", node_model="HoleGeneration", min_error=1e5) else: print("Avalanche probe skipped (AVALANCHE option is disabled).") print("--- END RECON PROBE ---\n") @@ -803,14 +837,48 @@ devsim.write_devices(file=f"{OUT_DIR}sweep_preview_final.tec", type="tecplot") np.savetxt(f"{OUT_DIR}sweep_iv_2d.csv", np.column_stack((voltage_list, current_list)), header="Voltage(V),Current(A/cm)", delimiter=",") +# Helper to read recon logs +def load_recon_log(filepath): + voltages = [] + currents = [] + if os.path.exists(filepath): + with open(filepath, "r") as f: + for line in f: + if "Voltage" in line or not line.strip(): + continue + parts = line.strip().split() + if len(parts) == 2 and parts[1] != "FAILED": + try: + v = float(parts[0]) + i = float(parts[1]) + voltages.append(v) + currents.append(abs(i)) + except ValueError: + pass + return voltages, currents + # Plot and save I-V curve plt.figure(figsize=(8, 6)) -plt.plot(voltage_list, np.abs(current_list), 'o-', color='#1f77b4', markersize=2) +plt.plot(voltage_list, np.abs(current_list), 'o-', color='#1f77b4', markersize=2, label='Main Sweep (Drift-Diffusion)') + +v_av, i_av = load_recon_log(f"{OUT_DIR}recon_avalanche.log") +if v_av: + plt.plot(v_av, i_av, 's--', color='#ff7f0e', markersize=2, label='Avalanche Only (Recon)') + +v_bt, i_bt = load_recon_log(f"{OUT_DIR}recon_btbt.log") +if v_bt: + plt.plot(v_bt, i_bt, '^--', color='#2ca02c', markersize=2, label='BTBT Only (Recon)') + +v_av_bt, i_av_bt = load_recon_log(f"{OUT_DIR}recon_av_btbt.log") +if v_av_bt: + plt.plot(v_av_bt, i_av_bt, 'd--', color='#d62728', markersize=2, label='Avalanche + BTBT (Recon)') + plt.yscale('log') plt.grid(True, which="both", ls="--") plt.xlabel("Bias Voltage (V)") plt.ylabel("Terminal Current (A/cm, Log Scale)") plt.title(f"{SIM_NAME} (T={temp_val}K)") +plt.legend() plt.tight_layout() plt.savefig(f"{OUT_DIR}sweep_iv_2d.png", dpi=300) plt.close() diff --git a/run_refinement_2d.py b/run_refinement_2d.py index 58571b8..a1407fe 100644 --- a/run_refinement_2d.py +++ b/run_refinement_2d.py @@ -2,8 +2,10 @@ import devsim import numpy as np import matplotlib.pyplot as plt import os +import os import sys -sys.path.append("/home/pchan/devsim2026") +DEV_DIR = os.environ.get("DEV_DIR", "devices/Triac_rp") +sys.path.insert(0, os.path.abspath(DEV_DIR)) from device_config import * from physics.model_create import * from physics.new_physics import * @@ -39,41 +41,45 @@ def run_simulation(mesh_file="device_2d.msh", tec_file="static_preview.tec", png devsim.create_device(mesh=device, device=device) # 2. Set up doping in Silicon region - devsim.node_model(device=device, region="Silicon", name="nD_sub", equation=f"{N_SUB}") - - def get_erfc_expr(peak, x1, x2, hdiff, vdiff): - return f"{peak} * erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))" - - p11_left_expr = get_erfc_expr(P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) - p11_right_expr = get_erfc_expr(P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) - devsim.node_model(device=device, region="Silicon", name="nA_p11_l", equation=p11_left_expr) - devsim.node_model(device=device, region="Silicon", name="nA_p11_r", equation=p11_right_expr) - - p12_left_expr = get_erfc_expr(P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) - p12_right_expr = get_erfc_expr(P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) - devsim.node_model(device=device, region="Silicon", name="nA_p12_l", equation=p12_left_expr) - devsim.node_model(device=device, region="Silicon", name="nA_p12_r", equation=p12_right_expr) - - p13_left_expr = get_erfc_expr(P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) - p13_right_expr = get_erfc_expr(P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) - devsim.node_model(device=device, region="Silicon", name="nA_p13_l", equation=p13_left_expr) - devsim.node_model(device=device, region="Silicon", name="nA_p13_r", equation=p13_right_expr) - - nplus_left_expr = get_erfc_expr(NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) - nplus_right_expr = get_erfc_expr(NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF) - devsim.node_model(device=device, region="Silicon", name="nD_nplus_l", equation=nplus_left_expr) - devsim.node_model(device=device, region="Silicon", name="nD_nplus_r", equation=nplus_right_expr) - - mring_l_expr = get_erfc_expr(NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) - mring_r_expr = get_erfc_expr(NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF) - devsim.node_model(device=device, region="Silicon", name="nD_mring_l", equation=mring_l_expr) - devsim.node_model(device=device, region="Silicon", name="nD_mring_r", equation=mring_r_expr) - - devsim.node_model(device=device, region="Silicon", name="Donors", - equation="nD_sub + nD_nplus_l + nD_nplus_r + nD_mring_l + nD_mring_r") - devsim.node_model(device=device, region="Silicon", name="Acceptors", - equation="1e10 + nA_p11_l + nA_p11_r + nA_p12_l + nA_p12_r + nA_p13_l + nA_p13_r") - devsim.node_model(device=device, region="Silicon", name="NetDoping", equation="Donors - Acceptors") + if os.environ.get("USE_PCAD", "false").lower() == "true": + from device_pcad_config import apply_pcad_doping_2d + apply_pcad_doping_2d(device, region="Silicon") + else: + devsim.node_model(device=device, region="Silicon", name="nD_sub", equation=f"{N_SUB}") + + def get_erfc_expr(peak, x1, x2, hdiff, vdiff): + return f"{peak} * erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))" + + p11_left_expr = get_erfc_expr(P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p11_right_expr = get_erfc_expr(P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p11_l", equation=p11_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p11_r", equation=p11_right_expr) + + p12_left_expr = get_erfc_expr(P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p12_right_expr = get_erfc_expr(P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p12_l", equation=p12_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p12_r", equation=p12_right_expr) + + p13_left_expr = get_erfc_expr(P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p13_right_expr = get_erfc_expr(P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p13_l", equation=p13_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p13_r", equation=p13_right_expr) + + nplus_left_expr = get_erfc_expr(NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) + nplus_right_expr = get_erfc_expr(NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nD_nplus_l", equation=nplus_left_expr) + devsim.node_model(device=device, region="Silicon", name="nD_nplus_r", equation=nplus_right_expr) + + mring_l_expr = get_erfc_expr(NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) + mring_r_expr = get_erfc_expr(NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nD_mring_l", equation=mring_l_expr) + devsim.node_model(device=device, region="Silicon", name="nD_mring_r", equation=mring_r_expr) + + devsim.node_model(device=device, region="Silicon", name="Donors", + equation="nD_sub + nD_nplus_l + nD_nplus_r + nD_mring_l + nD_mring_r") + devsim.node_model(device=device, region="Silicon", name="Acceptors", + equation="1e10 + nA_p11_l + nA_p11_r + nA_p12_l + nA_p12_r + nA_p13_l + nA_p13_r") + devsim.node_model(device=device, region="Silicon", name="NetDoping", equation="Donors - Acceptors") # 3. Solutions and Physics CreateSolution(device, "Silicon", "Potential") @@ -229,7 +235,10 @@ def run_simulation(mesh_file="device_2d.msh", tec_file="static_preview.tec", png def generate_background_mesh(): # 1. Run simulation on current mesh to get Emag - device = run_simulation("device_2d.msh", "static_preview.tec", "static_potential_2d.png", suffix="(Coarse Mesh)") + mesh_file = os.path.join(DEV_DIR, "device_2d.msh") + tec_file = os.path.join(DEV_DIR, "static_preview.tec") + png_file = os.path.join(DEV_DIR, "static_potential_2d.png") + device = run_simulation(mesh_file, tec_file, png_file, suffix="(Coarse Mesh)") # 2. Extract elements and Emag print("Generating background mesh...") @@ -240,7 +249,8 @@ def generate_background_mesh(): alpha = 1.0e-3 # Scaling coefficient for Emag # We will write to device_bgmesh.pos - with open("device_bgmesh.pos", "w") as f: + bgmesh_file = os.path.join(DEV_DIR, "device_bgmesh.pos") + with open(bgmesh_file, "w") as f: f.write('View "background mesh" {\n') # Write for Silicon, Oxide, Molding regions diff --git a/solve_static_2d.py b/solve_static_2d.py index 49f98eb..a95a190 100644 --- a/solve_static_2d.py +++ b/solve_static_2d.py @@ -10,7 +10,10 @@ os.environ["OPENBLAS_NUM_THREADS"] = "4" import devsim import numpy as np -OUT_DIR = "output_this_run/" +DEV_DIR = os.environ.get("DEV_DIR", "devices/Triac_rp") +sys.path.insert(0, os.path.abspath(DEV_DIR)) + +OUT_DIR = os.path.join(os.environ.get("OUT_DIR", os.path.join(DEV_DIR, "output_this_run")), "") os.makedirs(OUT_DIR, exist_ok=True) import matplotlib.pyplot as plt from device_config import * @@ -20,7 +23,8 @@ from physics.new_physics import * device = "device_2d" # 1. Load the mesh -devsim.create_gmsh_mesh(mesh=device, file="device_2d.msh") +mesh_file = os.path.join(DEV_DIR, "device_2d.msh") +devsim.create_gmsh_mesh(mesh=device, file=mesh_file) devsim.add_gmsh_region(mesh=device, gmsh_name="Silicon", region="Silicon", material="Silicon") devsim.add_gmsh_region(mesh=device, gmsh_name="Oxide", region="Oxide", material="Oxide") devsim.add_gmsh_region(mesh=device, gmsh_name="Molding", region="Molding", material="Molding") @@ -51,46 +55,50 @@ devsim.create_device(mesh=device, device=device) # 2. Set up doping in Silicon region -devsim.node_model(device=device, region="Silicon", name="nD_sub", equation=f"{N_SUB}") +if os.environ.get("USE_PCAD", "false").lower() == "true": + from device_pcad_config import apply_pcad_doping_2d + apply_pcad_doping_2d(device, region="Silicon") +else: + devsim.node_model(device=device, region="Silicon", name="nD_sub", equation=f"{N_SUB}") -def get_erfc_expr(peak, x1, x2, hdiff, vdiff): - return f"{peak} * erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))" + def get_erfc_expr(peak, x1, x2, hdiff, vdiff): + return f"{peak} * erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))" -# P-wells -p11_left_expr = get_erfc_expr(P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p11_right_expr = get_erfc_expr(P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p11_l", equation=p11_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p11_r", equation=p11_right_expr) + # P-wells + p11_left_expr = get_erfc_expr(P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p11_right_expr = get_erfc_expr(P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p11_l", equation=p11_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p11_r", equation=p11_right_expr) -p12_left_expr = get_erfc_expr(P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p12_right_expr = get_erfc_expr(P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p12_l", equation=p12_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p12_r", equation=p12_right_expr) + p12_left_expr = get_erfc_expr(P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p12_right_expr = get_erfc_expr(P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p12_l", equation=p12_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p12_r", equation=p12_right_expr) -p13_left_expr = get_erfc_expr(P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p13_right_expr = get_erfc_expr(P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p13_l", equation=p13_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p13_r", equation=p13_right_expr) + p13_left_expr = get_erfc_expr(P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p13_right_expr = get_erfc_expr(P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p13_l", equation=p13_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p13_r", equation=p13_right_expr) -# N+ -nplus_left_expr = get_erfc_expr(NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) -nplus_right_expr = get_erfc_expr(NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nD_nplus_l", equation=nplus_left_expr) -devsim.node_model(device=device, region="Silicon", name="nD_nplus_r", equation=nplus_right_expr) + # N+ + nplus_left_expr = get_erfc_expr(NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) + nplus_right_expr = get_erfc_expr(NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nD_nplus_l", equation=nplus_left_expr) + devsim.node_model(device=device, region="Silicon", name="nD_nplus_r", equation=nplus_right_expr) -# MRING -mring_l_expr = get_erfc_expr(NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) -mring_r_expr = get_erfc_expr(NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nD_mring_l", equation=mring_l_expr) -devsim.node_model(device=device, region="Silicon", name="nD_mring_r", equation=mring_r_expr) + # MRING + mring_l_expr = get_erfc_expr(NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) + mring_r_expr = get_erfc_expr(NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nD_mring_l", equation=mring_l_expr) + devsim.node_model(device=device, region="Silicon", name="nD_mring_r", equation=mring_r_expr) -# Combine into Donors and Acceptors -devsim.node_model(device=device, region="Silicon", name="Donors", - equation="nD_sub + nD_nplus_l + nD_nplus_r + nD_mring_l + nD_mring_r") -devsim.node_model(device=device, region="Silicon", name="Acceptors", - equation="1e10 + nA_p11_l + nA_p11_r + nA_p12_l + nA_p12_r + nA_p13_l + nA_p13_r") -devsim.node_model(device=device, region="Silicon", name="NetDoping", equation="Donors - Acceptors") -devsim.node_model(device=device, region="Silicon", name="LogNetDoping", equation="asinh(NetDoping / 2.0) / log(10.0)") + # Combine into Donors and Acceptors + devsim.node_model(device=device, region="Silicon", name="Donors", + equation="nD_sub + nD_nplus_l + nD_nplus_r + nD_mring_l + nD_mring_r") + devsim.node_model(device=device, region="Silicon", name="Acceptors", + equation="1e10 + nA_p11_l + nA_p11_r + nA_p12_l + nA_p12_r + nA_p13_l + nA_p13_r") + devsim.node_model(device=device, region="Silicon", name="NetDoping", equation="Donors - Acceptors") + devsim.node_model(device=device, region="Silicon", name="LogNetDoping", equation="asinh(NetDoping / 2.0) / log(10.0)") # 3. Create solution variables and physics models CreateSolution(device, "Silicon", "Potential") diff --git a/solve_sweep_2d.py b/solve_sweep_2d.py index 40e4e92..5f8e6c9 100644 --- a/solve_sweep_2d.py +++ b/solve_sweep_2d.py @@ -19,7 +19,10 @@ if mkl_libs: import devsim import numpy as np -OUT_DIR = "output_this_run/" +DEV_DIR = os.environ.get("DEV_DIR", "devices/Triac_rp") +sys.path.insert(0, os.path.abspath(DEV_DIR)) + +OUT_DIR = os.path.join(os.environ.get("OUT_DIR", os.path.join(DEV_DIR, "output_this_run")), "") os.makedirs(OUT_DIR, exist_ok=True) import matplotlib.pyplot as plt import time @@ -27,7 +30,6 @@ import time # Enable Intel MKL PARDISO multi-threaded sparse solver devsim.set_parameter(name="solver_type", value="pardiso") -sys.path.append("/home/pchan/devsim2026") from device_config import * from physics.model_create import * from physics.new_physics import * @@ -35,8 +37,9 @@ from physics.new_physics import * device = "device_2d" # 1. Load the mesh -print("Loading mesh: device_2d.msh...") -devsim.create_gmsh_mesh(mesh=device, file="device_2d.msh") +mesh_file = os.path.join(DEV_DIR, "device_2d.msh") +print(f"Loading mesh: {mesh_file}...") +devsim.create_gmsh_mesh(mesh=device, file=mesh_file) devsim.add_gmsh_region(mesh=device, gmsh_name="Silicon", region="Silicon", material="Silicon") devsim.add_gmsh_region(mesh=device, gmsh_name="Oxide", region="Oxide", material="Oxide") devsim.add_gmsh_region(mesh=device, gmsh_name="Molding", region="Molding", material="Molding") @@ -64,42 +67,46 @@ devsim.finalize_mesh(mesh=device) devsim.create_device(mesh=device, device=device) # 2. Set up doping in Silicon region -devsim.node_model(device=device, region="Silicon", name="nD_sub", equation=f"{N_SUB}") +if os.environ.get("USE_PCAD", "false").lower() == "true": + from device_pcad_config import apply_pcad_doping_2d + apply_pcad_doping_2d(device, region="Silicon") +else: + devsim.node_model(device=device, region="Silicon", name="nD_sub", equation=f"{N_SUB}") -def get_erfc_expr(peak, x1, x2, hdiff, vdiff): - return f"{peak} * erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))" + def get_erfc_expr(peak, x1, x2, hdiff, vdiff): + return f"{peak} * erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))" -p11_left_expr = get_erfc_expr(P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p11_right_expr = get_erfc_expr(P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p11_l", equation=p11_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p11_r", equation=p11_right_expr) + p11_left_expr = get_erfc_expr(P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p11_right_expr = get_erfc_expr(P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p11_l", equation=p11_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p11_r", equation=p11_right_expr) -p12_left_expr = get_erfc_expr(P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p12_right_expr = get_erfc_expr(P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p12_l", equation=p12_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p12_r", equation=p12_right_expr) + p12_left_expr = get_erfc_expr(P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p12_right_expr = get_erfc_expr(P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p12_l", equation=p12_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p12_r", equation=p12_right_expr) -p13_left_expr = get_erfc_expr(P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p13_right_expr = get_erfc_expr(P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p13_l", equation=p13_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p13_r", equation=p13_right_expr) + p13_left_expr = get_erfc_expr(P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p13_right_expr = get_erfc_expr(P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p13_l", equation=p13_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p13_r", equation=p13_right_expr) -nplus_left_expr = get_erfc_expr(NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) -nplus_right_expr = get_erfc_expr(NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nD_nplus_l", equation=nplus_left_expr) -devsim.node_model(device=device, region="Silicon", name="nD_nplus_r", equation=nplus_right_expr) + nplus_left_expr = get_erfc_expr(NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) + nplus_right_expr = get_erfc_expr(NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nD_nplus_l", equation=nplus_left_expr) + devsim.node_model(device=device, region="Silicon", name="nD_nplus_r", equation=nplus_right_expr) -mring_l_expr = get_erfc_expr(NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) -mring_r_expr = get_erfc_expr(NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nD_mring_l", equation=mring_l_expr) -devsim.node_model(device=device, region="Silicon", name="nD_mring_r", equation=mring_r_expr) + mring_l_expr = get_erfc_expr(NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) + mring_r_expr = get_erfc_expr(NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nD_mring_l", equation=mring_l_expr) + devsim.node_model(device=device, region="Silicon", name="nD_mring_r", equation=mring_r_expr) -devsim.node_model(device=device, region="Silicon", name="Donors", - equation="nD_sub + nD_nplus_l + nD_nplus_r + nD_mring_l + nD_mring_r") -devsim.node_model(device=device, region="Silicon", name="Acceptors", - equation="1e10 + nA_p11_l + nA_p11_r + nA_p12_l + nA_p12_r + nA_p13_l + nA_p13_r") -devsim.node_model(device=device, region="Silicon", name="NetDoping", equation="Donors - Acceptors") -devsim.node_model(device=device, region="Silicon", name="LogNetDoping", equation="asinh(NetDoping / 2.0) / log(10.0)") + devsim.node_model(device=device, region="Silicon", name="Donors", + equation="nD_sub + nD_nplus_l + nD_nplus_r + nD_mring_l + nD_mring_r") + devsim.node_model(device=device, region="Silicon", name="Acceptors", + equation="1e10 + nA_p11_l + nA_p11_r + nA_p12_l + nA_p12_r + nA_p13_l + nA_p13_r") + devsim.node_model(device=device, region="Silicon", name="NetDoping", equation="Donors - Acceptors") + devsim.node_model(device=device, region="Silicon", name="LogNetDoping", equation="asinh(NetDoping / 2.0) / log(10.0)") # 3. Initialize electrostatic potential simulation (Poisson only) CreateSolution(device, "Silicon", "Potential") diff --git a/solve_sweep_bv.py b/solve_sweep_bv.py index 88911b8..8cfd264 100644 --- a/solve_sweep_bv.py +++ b/solve_sweep_bv.py @@ -21,7 +21,10 @@ if mkl_libs: import devsim import numpy as np -OUT_DIR = "output_this_run/" +DEV_DIR = os.environ.get("DEV_DIR", "devices/Triac_rp") +sys.path.insert(0, os.path.abspath(DEV_DIR)) + +OUT_DIR = os.path.join(os.environ.get("OUT_DIR", os.path.join(DEV_DIR, "output_this_run")), "") os.makedirs(OUT_DIR, exist_ok=True) import matplotlib.pyplot as plt import time @@ -29,7 +32,6 @@ import time # Enable Intel MKL PARDISO multi-threaded sparse solver devsim.set_parameter(name="solver_type", value="pardiso") -sys.path.append("/home/pchan/devsim2026") from device_config import * from physics.model_create import * from physics.new_physics import * @@ -37,8 +39,9 @@ from physics.new_physics import * device = "device_2d" # 1. Load the mesh -print("Loading mesh: device_2d.msh...") -devsim.create_gmsh_mesh(mesh=device, file="device_2d.msh") +mesh_file = os.path.join(DEV_DIR, "device_2d.msh") +print(f"Loading mesh: {mesh_file}...") +devsim.create_gmsh_mesh(mesh=device, file=mesh_file) devsim.add_gmsh_region(mesh=device, gmsh_name="Silicon", region="Silicon", material="Silicon") devsim.add_gmsh_region(mesh=device, gmsh_name="Oxide", region="Oxide", material="Oxide") devsim.add_gmsh_region(mesh=device, gmsh_name="Molding", region="Molding", material="Molding") @@ -66,42 +69,46 @@ devsim.finalize_mesh(mesh=device) devsim.create_device(mesh=device, device=device) # 2. Set up doping in Silicon region -devsim.node_model(device=device, region="Silicon", name="nD_sub", equation=f"{N_SUB}") +if os.environ.get("USE_PCAD", "false").lower() == "true": + from device_pcad_config import apply_pcad_doping_2d + apply_pcad_doping_2d(device, region="Silicon") +else: + devsim.node_model(device=device, region="Silicon", name="nD_sub", equation=f"{N_SUB}") -def get_erfc_expr(peak, x1, x2, hdiff, vdiff): - return f"{peak} * erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))" + def get_erfc_expr(peak, x1, x2, hdiff, vdiff): + return f"{peak} * erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))" -p11_left_expr = get_erfc_expr(P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p11_right_expr = get_erfc_expr(P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p11_l", equation=p11_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p11_r", equation=p11_right_expr) + p11_left_expr = get_erfc_expr(P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p11_right_expr = get_erfc_expr(P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p11_l", equation=p11_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p11_r", equation=p11_right_expr) -p12_left_expr = get_erfc_expr(P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p12_right_expr = get_erfc_expr(P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p12_l", equation=p12_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p12_r", equation=p12_right_expr) + p12_left_expr = get_erfc_expr(P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p12_right_expr = get_erfc_expr(P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p12_l", equation=p12_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p12_r", equation=p12_right_expr) -p13_left_expr = get_erfc_expr(P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p13_right_expr = get_erfc_expr(P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p13_l", equation=p13_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p13_r", equation=p13_right_expr) + p13_left_expr = get_erfc_expr(P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p13_right_expr = get_erfc_expr(P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p13_l", equation=p13_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p13_r", equation=p13_right_expr) -nplus_left_expr = get_erfc_expr(NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) -nplus_right_expr = get_erfc_expr(NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nD_nplus_l", equation=nplus_left_expr) -devsim.node_model(device=device, region="Silicon", name="nD_nplus_r", equation=nplus_right_expr) + nplus_left_expr = get_erfc_expr(NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) + nplus_right_expr = get_erfc_expr(NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nD_nplus_l", equation=nplus_left_expr) + devsim.node_model(device=device, region="Silicon", name="nD_nplus_r", equation=nplus_right_expr) -mring_l_expr = get_erfc_expr(NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) -mring_r_expr = get_erfc_expr(NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nD_mring_l", equation=mring_l_expr) -devsim.node_model(device=device, region="Silicon", name="nD_mring_r", equation=mring_r_expr) + mring_l_expr = get_erfc_expr(NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) + mring_r_expr = get_erfc_expr(NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nD_mring_l", equation=mring_l_expr) + devsim.node_model(device=device, region="Silicon", name="nD_mring_r", equation=mring_r_expr) -devsim.node_model(device=device, region="Silicon", name="Donors", - equation="nD_sub + nD_nplus_l + nD_nplus_r + nD_mring_l + nD_mring_r") -devsim.node_model(device=device, region="Silicon", name="Acceptors", - equation="1e10 + nA_p11_l + nA_p11_r + nA_p12_l + nA_p12_r + nA_p13_l + nA_p13_r") -devsim.node_model(device=device, region="Silicon", name="NetDoping", equation="Donors - Acceptors") -devsim.node_model(device=device, region="Silicon", name="LogNetDoping", equation="asinh(NetDoping / 2.0) / log(10.0)") + devsim.node_model(device=device, region="Silicon", name="Donors", + equation="nD_sub + nD_nplus_l + nD_nplus_r + nD_mring_l + nD_mring_r") + devsim.node_model(device=device, region="Silicon", name="Acceptors", + equation="1e10 + nA_p11_l + nA_p11_r + nA_p12_l + nA_p12_r + nA_p13_l + nA_p13_r") + devsim.node_model(device=device, region="Silicon", name="NetDoping", equation="Donors - Acceptors") + devsim.node_model(device=device, region="Silicon", name="LogNetDoping", equation="asinh(NetDoping / 2.0) / log(10.0)") # 3. Initialize electrostatic potential simulation (Poisson only) CreateSolution(device, "Silicon", "Potential") diff --git a/solve_sweep_recon.py b/solve_sweep_recon.py index f786f02..f5605c8 100644 --- a/solve_sweep_recon.py +++ b/solve_sweep_recon.py @@ -20,7 +20,10 @@ if mkl_libs: import devsim import numpy as np -OUT_DIR = "output_this_run/" +DEV_DIR = os.environ.get("DEV_DIR", "devices/Triac_rp") +sys.path.insert(0, os.path.abspath(DEV_DIR)) + +OUT_DIR = os.path.join(os.environ.get("OUT_DIR", os.path.join(DEV_DIR, "output_this_run")), "") os.makedirs(OUT_DIR, exist_ok=True) import matplotlib.pyplot as plt import time @@ -28,7 +31,6 @@ import time # Enable Intel MKL PARDISO multi-threaded sparse solver devsim.set_parameter(name="solver_type", value="pardiso") -sys.path.append("/home/pchan/devsim2026") from device_config import * from physics.model_create import * from physics.new_physics import * @@ -36,8 +38,9 @@ from physics.new_physics import * device = "device_2d" # 1. Load the mesh -print("Loading mesh: device_2d.msh...") -devsim.create_gmsh_mesh(mesh=device, file="device_2d.msh") +mesh_file = os.path.join(DEV_DIR, "device_2d.msh") +print(f"Loading mesh: {mesh_file}...") +devsim.create_gmsh_mesh(mesh=device, file=mesh_file) devsim.add_gmsh_region(mesh=device, gmsh_name="Silicon", region="Silicon", material="Silicon") devsim.add_gmsh_region(mesh=device, gmsh_name="Oxide", region="Oxide", material="Oxide") devsim.add_gmsh_region(mesh=device, gmsh_name="Molding", region="Molding", material="Molding") @@ -65,42 +68,46 @@ devsim.finalize_mesh(mesh=device) devsim.create_device(mesh=device, device=device) # 2. Set up doping in Silicon region -devsim.node_model(device=device, region="Silicon", name="nD_sub", equation=f"{N_SUB}") +if os.environ.get("USE_PCAD", "false").lower() == "true": + from device_pcad_config import apply_pcad_doping_2d + apply_pcad_doping_2d(device, region="Silicon") +else: + devsim.node_model(device=device, region="Silicon", name="nD_sub", equation=f"{N_SUB}") -def get_erfc_expr(peak, x1, x2, hdiff, vdiff): - return f"{peak} * erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))" + def get_erfc_expr(peak, x1, x2, hdiff, vdiff): + return f"{peak} * erfc(y / {vdiff}) * 0.5 * (erf((x - ({x1})) / {hdiff}) - erf((x - ({x2})) / {hdiff}))" -p11_left_expr = get_erfc_expr(P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p11_right_expr = get_erfc_expr(P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p11_l", equation=p11_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p11_r", equation=p11_right_expr) + p11_left_expr = get_erfc_expr(P11_PEAK, -P11_X2, -P11_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p11_right_expr = get_erfc_expr(P11_PEAK, P11_X1, P11_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p11_l", equation=p11_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p11_r", equation=p11_right_expr) -p12_left_expr = get_erfc_expr(P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p12_right_expr = get_erfc_expr(P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p12_l", equation=p12_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p12_r", equation=p12_right_expr) + p12_left_expr = get_erfc_expr(P12_PEAK, -P12_X2, -P12_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p12_right_expr = get_erfc_expr(P12_PEAK, P12_X1, P12_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p12_l", equation=p12_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p12_r", equation=p12_right_expr) -p13_left_expr = get_erfc_expr(P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) -p13_right_expr = get_erfc_expr(P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nA_p13_l", equation=p13_left_expr) -devsim.node_model(device=device, region="Silicon", name="nA_p13_r", equation=p13_right_expr) + p13_left_expr = get_erfc_expr(P13_PEAK, -P13_X2, -P13_X1, P_WELL_HDDIFF, P_WELL_VDDIFF) + p13_right_expr = get_erfc_expr(P13_PEAK, P13_X1, P13_X2, P_WELL_HDDIFF, P_WELL_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nA_p13_l", equation=p13_left_expr) + devsim.node_model(device=device, region="Silicon", name="nA_p13_r", equation=p13_right_expr) -nplus_left_expr = get_erfc_expr(NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) -nplus_right_expr = get_erfc_expr(NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nD_nplus_l", equation=nplus_left_expr) -devsim.node_model(device=device, region="Silicon", name="nD_nplus_r", equation=nplus_right_expr) + nplus_left_expr = get_erfc_expr(NPLUS_PEAK, -NPLUS_X2, -NPLUS_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) + nplus_right_expr = get_erfc_expr(NPLUS_PEAK, NPLUS_X1, NPLUS_X2, NPLUS_HDDIFF, NPLUS_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nD_nplus_l", equation=nplus_left_expr) + devsim.node_model(device=device, region="Silicon", name="nD_nplus_r", equation=nplus_right_expr) -mring_l_expr = get_erfc_expr(NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) -mring_r_expr = get_erfc_expr(NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF) -devsim.node_model(device=device, region="Silicon", name="nD_mring_l", equation=mring_l_expr) -devsim.node_model(device=device, region="Silicon", name="nD_mring_r", equation=mring_r_expr) + mring_l_expr = get_erfc_expr(NPLUS_PEAK, -W_DEVICE, -MRING_X1, NPLUS_HDDIFF, NPLUS_VDDIFF) + mring_r_expr = get_erfc_expr(NPLUS_PEAK, MRING_X1, W_DEVICE, NPLUS_HDDIFF, NPLUS_VDDIFF) + devsim.node_model(device=device, region="Silicon", name="nD_mring_l", equation=mring_l_expr) + devsim.node_model(device=device, region="Silicon", name="nD_mring_r", equation=mring_r_expr) -devsim.node_model(device=device, region="Silicon", name="Donors", - equation="nD_sub + nD_nplus_l + nD_nplus_r + nD_mring_l + nD_mring_r") -devsim.node_model(device=device, region="Silicon", name="Acceptors", - equation="1e10 + nA_p11_l + nA_p11_r + nA_p12_l + nA_p12_r + nA_p13_l + nA_p13_r") -devsim.node_model(device=device, region="Silicon", name="NetDoping", equation="Donors - Acceptors") -devsim.node_model(device=device, region="Silicon", name="LogNetDoping", equation="asinh(NetDoping / 2.0) / log(10.0)") + devsim.node_model(device=device, region="Silicon", name="Donors", + equation="nD_sub + nD_nplus_l + nD_nplus_r + nD_mring_l + nD_mring_r") + devsim.node_model(device=device, region="Silicon", name="Acceptors", + equation="1e10 + nA_p11_l + nA_p11_r + nA_p12_l + nA_p12_r + nA_p13_l + nA_p13_r") + devsim.node_model(device=device, region="Silicon", name="NetDoping", equation="Donors - Acceptors") + devsim.node_model(device=device, region="Silicon", name="LogNetDoping", equation="asinh(NetDoping / 2.0) / log(10.0)") # 3. Initialize electrostatic potential simulation (Poisson only) CreateSolution(device, "Silicon", "Potential") @@ -290,7 +297,7 @@ current_list = [0.0] # Recon variables next_recon_v = 50.0 -with open("recon_avalanche.log", "w") as f: +with open(f"{OUT_DIR}recon_avalanche.log", "w") as f: f.write("Voltage(V)\tAvalancheCurrent(A)\n") @@ -383,7 +390,7 @@ while v_current < v_target: if v_current >= next_recon_v: state_data = save_state(device) seed_data = {"voltage": v_current, "step_size": step_size, "state": state_data} - seed_filename = f"seed_{int(next_recon_v)}V.pkl" + seed_filename = f"{OUT_DIR}seed_{int(next_recon_v)}V.pkl" with open(seed_filename, "wb") as f: pickle.dump(seed_data, f) print(f"\n--- RECON PROBE at {v_current:.2f} V ---") @@ -413,15 +420,15 @@ while v_current < v_target: ia_p_p12 = devsim.get_contact_current(device=device, contact="MT1_P12_Si", equation="HoleContinuityEquation") av_curr = ia_n_si + ia_p_si + ia_n_p12 + ia_p_p12 print(f"Avalanche Current at {v_current:.2f} V: {av_curr:.4e} A") - with open("recon_avalanche.log", "a") as f: + with open(f"{OUT_DIR}recon_avalanche.log", "a") as f: f.write(f"{v_current:.2f}\t{av_curr:.4e}\n") else: print("Avalanche failed to converge.") - with open("recon_avalanche.log", "a") as f: + with open(f"{OUT_DIR}recon_avalanche.log", "a") as f: f.write(f"{v_current:.2f}\tFAILED\n") except devsim.error: print("Avalanche failed to converge.") - with open("recon_avalanche.log", "a") as f: + with open(f"{OUT_DIR}recon_avalanche.log", "a") as f: f.write(f"{v_current:.2f}\tFAILED\n") # Restore state and Turn OFF Avalanche