pchang718
371 lines
18 KiB
Python
371 lines
18 KiB
Python
# mos_transfer_single_temp.py
|
|
# 2D MOSFET Ids-Vgs Transfer Curve Sweep at a Specific Temperature
|
|
# Sweeps Vgs from 0.0V to 10.0V under constant Vds = 1.0V
|
|
|
|
import os
|
|
import sys
|
|
import glob
|
|
import gc
|
|
import time
|
|
import argparse
|
|
import pickle
|
|
import numpy as np
|
|
|
|
# Limit thread count to avoid system resource starvation
|
|
os.environ["OMP_NUM_THREADS"] = "4"
|
|
os.environ["MKL_NUM_THREADS"] = "4"
|
|
os.environ["TBB_NUM_THREADS"] = "4"
|
|
os.environ["OPENBLAS_NUM_THREADS"] = "4"
|
|
os.environ["MKL_DISABLE_FAST_MM"] = "1"
|
|
|
|
# Load Intel MKL runtime library if available
|
|
mkl_libs = glob.glob(os.path.join(os.path.dirname(sys.executable), "../lib/libmkl_rt.so.*"))
|
|
if mkl_libs:
|
|
os.environ["DEVSIM_MATH_LIBS"] = mkl_libs[0]
|
|
|
|
import devsim
|
|
|
|
# Use MKL PARDISO solver
|
|
devsim.set_parameter(name="solver_type", value="pardiso")
|
|
|
|
# Parse arguments
|
|
parser = argparse.ArgumentParser(description="Run transfer characteristics sweep at a single temperature.")
|
|
parser.add_argument("--temp", type=float, required=True, help="Temperature in Celsius.")
|
|
parser.add_argument("--out_dir", type=str, required=True, help="Output directory path.")
|
|
parser.add_argument("--vds", type=float, default=1.0, help="Drain-Source Voltage Vds.")
|
|
args = parser.parse_args()
|
|
|
|
temp_C = args.temp
|
|
T_kelvin = temp_C + 273.15
|
|
OUT_DIR = os.path.join(args.out_dir, "")
|
|
os.makedirs(OUT_DIR, exist_ok=True)
|
|
|
|
# Default device directory to LDMOS
|
|
DEV_DIR = os.environ.get("DEV_DIR", "devices/LDMOS")
|
|
sys.path.insert(0, os.path.abspath(DEV_DIR))
|
|
|
|
# Import geometry parameters and physics creators
|
|
from device_config import *
|
|
from physics.model_create import *
|
|
from physics.new_physics import *
|
|
|
|
device = "device_2d"
|
|
|
|
# --- 1. Load Mesh ---
|
|
mesh_file = os.path.join(DEV_DIR, "device_2d.msh")
|
|
print(f"[{temp_C} C] 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")
|
|
|
|
# Add contacts
|
|
devsim.add_gmsh_contact(mesh=device, gmsh_name="contR_Si", name="contR_Si", region="Silicon", material="metal")
|
|
devsim.add_gmsh_contact(mesh=device, gmsh_name="contL_Si", name="contL_Si", region="Silicon", material="metal")
|
|
devsim.add_gmsh_contact(mesh=device, gmsh_name="gate_Ox", name="gate_Ox", region="Oxide", material="metal")
|
|
|
|
devsim.add_gmsh_contact(mesh=device, gmsh_name="gate_Mold", name="gate_Mold", region="Molding", material="metal")
|
|
devsim.add_gmsh_contact(mesh=device, gmsh_name="contR_Mold", name="contR_Mold", region="Molding", material="metal")
|
|
devsim.add_gmsh_contact(mesh=device, gmsh_name="contL_Mold", name="contL_Mold", region="Molding", material="metal")
|
|
|
|
# Interfaces
|
|
devsim.add_gmsh_interface(mesh=device, gmsh_name="Si_Ox_Interface", name="Si_Ox", region0="Silicon", region1="Oxide")
|
|
devsim.add_gmsh_interface(mesh=device, gmsh_name="Ox_Mold_Interface", name="Ox_Mold", region0="Oxide", region1="Molding")
|
|
|
|
devsim.finalize_mesh(mesh=device)
|
|
devsim.create_device(mesh=device, device=device)
|
|
|
|
# --- 2. Setup Doping (PCAD) ---
|
|
try:
|
|
from device_pcad_config import apply_pcad_doping_2d
|
|
apply_pcad_doping_2d(device, region="Silicon")
|
|
print(f"[{temp_C} C] Applied PCAD doping successfully.")
|
|
except ImportError:
|
|
print(f"[{temp_C} C] Error: device_pcad_config not found. Make sure DEV_DIR points to the correct LDMOS folder.")
|
|
sys.exit(1)
|
|
|
|
# --- 3. Initialize Electrostatic Solution (Poisson) ---
|
|
CreateSolution(device, "Silicon", "Potential")
|
|
CreateSiliconPotentialOnly(device, "Silicon")
|
|
|
|
# Set temperature T on both device level and region level to prevent SetSiliconParameters overriding it back to 300K
|
|
devsim.set_parameter(device=device, name="T", value=T_kelvin)
|
|
devsim.set_parameter(device=device, region="Silicon", name="T", value=T_kelvin)
|
|
|
|
|
|
def CreateOxidePotentialOnly(device, region):
|
|
if not InNodeModelList(device, region, "Potential"):
|
|
CreateSolution(device, region, "Potential")
|
|
devsim.set_parameter(device=device, region=region, name="Permittivity", value=3.9 * 8.85e-14)
|
|
efield = "(Potential@n0 - Potential@n1)*EdgeInverseLength"
|
|
CreateEdgeModel(device, region, "EField", efield)
|
|
CreateEdgeModelDerivatives(device, region, "EField", efield, "Potential")
|
|
dfield = "Permittivity*EField"
|
|
CreateEdgeModel(device, region, "PotentialEdgeFlux", dfield)
|
|
CreateEdgeModelDerivatives(device, region, "PotentialEdgeFlux", dfield, "Potential")
|
|
devsim.equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
|
|
edge_model="PotentialEdgeFlux", variable_update="default", min_error=1e-3)
|
|
|
|
CreateOxidePotentialOnly(device, "Oxide")
|
|
|
|
def CreateMoldingPotentialOnly(device, region):
|
|
if not InNodeModelList(device, region, "Potential"):
|
|
CreateSolution(device, region, "Potential")
|
|
devsim.set_parameter(device=device, region=region, name="Permittivity", value=4.0 * 8.85e-14)
|
|
efield = "(Potential@n0 - Potential@n1)*EdgeInverseLength"
|
|
CreateEdgeModel(device, region, "EField", efield)
|
|
CreateEdgeModelDerivatives(device, region, "EField", efield, "Potential")
|
|
dfield = "Permittivity*EField"
|
|
CreateEdgeModel(device, region, "PotentialEdgeFlux", dfield)
|
|
CreateEdgeModelDerivatives(device, region, "PotentialEdgeFlux", dfield, "Potential")
|
|
devsim.equation(device=device, region=region, name="PotentialEquation", variable_name="Potential",
|
|
edge_model="PotentialEdgeFlux", variable_update="default", min_error=1e-3)
|
|
|
|
CreateMoldingPotentialOnly(device, "Molding")
|
|
|
|
# Interfaces continuous potential
|
|
def CreateContinuousPotentialInterface(device, interface):
|
|
model_name = CreateContinuousInterfaceModel(device, interface, "Potential")
|
|
devsim.interface_equation(device=device, interface=interface, name="PotentialEquation",
|
|
interface_model=model_name, type="continuous")
|
|
|
|
CreateContinuousPotentialInterface(device, "Si_Ox")
|
|
CreateContinuousPotentialInterface(device, "Ox_Mold")
|
|
|
|
# Set contact potential boundary equations
|
|
silicon_contacts = ["contL_Si", "contR_Si"]
|
|
for c in silicon_contacts:
|
|
devsim.set_parameter(device=device, name=GetContactBiasName(c), value=0.0)
|
|
CreateSiliconPotentialOnlyContact(device, "Silicon", c)
|
|
|
|
def CreateOxidePotentialOnlyContact(device, region, contact):
|
|
contact_bias = GetContactBiasName(contact)
|
|
contact_model = f"Potential - {contact_bias}"
|
|
contact_model_name = f"{contact}_bc"
|
|
CreateContactNodeModel(device, contact, contact_model_name, contact_model)
|
|
CreateContactNodeModelDerivative(device, contact, contact_model_name, contact_model, "Potential")
|
|
devsim.contact_equation(device=device, contact=contact, name="PotentialEquation",
|
|
node_model=contact_model_name, edge_charge_model="PotentialEdgeFlux")
|
|
|
|
for c in ["gate_Ox"]:
|
|
devsim.set_parameter(device=device, name=GetContactBiasName(c), value=0.0)
|
|
CreateOxidePotentialOnlyContact(device, "Oxide", c)
|
|
|
|
def CreateMoldingPotentialOnlyContact(device, region, contact):
|
|
contact_bias = GetContactBiasName(contact)
|
|
contact_model = f"Potential - {contact_bias}"
|
|
contact_model_name = f"{contact}_bc"
|
|
CreateContactNodeModel(device, contact, contact_model_name, contact_model)
|
|
CreateContactNodeModelDerivative(device, contact, contact_model_name, contact_model, "Potential")
|
|
devsim.contact_equation(device=device, contact=contact, name="PotentialEquation",
|
|
node_model=contact_model_name, edge_charge_model="PotentialEdgeFlux")
|
|
|
|
for c in ["gate_Mold", "contR_Mold", "contL_Mold"]:
|
|
devsim.set_parameter(device=device, name=GetContactBiasName(c), value=0.0)
|
|
CreateMoldingPotentialOnlyContact(device, "Molding", c)
|
|
|
|
# Solve initial zero-bias Poisson
|
|
print(f"[{temp_C} C] Solving initial Poisson at thermal equilibrium...")
|
|
devsim.solve(type="dc", absolute_error=1.0, relative_error=1e-10, maximum_iterations=50)
|
|
print(f"[{temp_C} C] Initial Poisson converged.")
|
|
|
|
# --- 4. Setup Drift-Diffusion (Poisson + Continuity) ---
|
|
CreateSolution(device, "Silicon", "Electrons")
|
|
CreateSolution(device, "Silicon", "Holes")
|
|
devsim.set_node_values(device=device, region="Silicon", name="Electrons", init_from="IntrinsicElectrons")
|
|
devsim.set_node_values(device=device, region="Silicon", name="Holes", init_from="IntrinsicHoles")
|
|
|
|
print(f"[{temp_C} C] Redefining equilibrium models to prevent high-bias overflow...")
|
|
devsim.node_model(device=device, region="Silicon", name="IntrinsicElectrons", equation="Electrons")
|
|
devsim.node_model(device=device, region="Silicon", name="IntrinsicElectrons:Potential", equation="0")
|
|
devsim.node_model(device=device, region="Silicon", name="IntrinsicElectrons:Electrons", equation="1")
|
|
devsim.node_model(device=device, region="Silicon", name="IntrinsicElectrons:Holes", equation="0")
|
|
|
|
devsim.node_model(device=device, region="Silicon", name="IntrinsicHoles", equation="Holes")
|
|
devsim.node_model(device=device, region="Silicon", name="IntrinsicHoles:Potential", equation="0")
|
|
devsim.node_model(device=device, region="Silicon", name="IntrinsicHoles:Electrons", equation="0")
|
|
devsim.node_model(device=device, region="Silicon", name="IntrinsicHoles:Holes", equation="1")
|
|
|
|
devsim.node_model(device=device, region="Silicon", name="PotentialIntrinsicCharge", equation="0")
|
|
devsim.node_model(device=device, region="Silicon", name="PotentialIntrinsicCharge:Potential", equation="0")
|
|
|
|
opts = CreateAroraMobilityLF(device, "Silicon")
|
|
CreateSiliconDriftDiffusion(device, "Silicon", **opts)
|
|
CreateAvalancheGeneration(device, "Silicon", opts['Jn'], opts['Jp'])
|
|
|
|
for c in silicon_contacts:
|
|
CreateSiliconDriftDiffusionContact(device, "Silicon", c, opts['Jn'], opts['Jp'])
|
|
|
|
# Set up positive update continuity equations for stable sweeping
|
|
devsim.equation(device=device, region="Silicon", name="ElectronContinuityEquation", variable_name="Electrons",
|
|
time_node_model="NCharge", edge_model=opts['Jn'], edge_volume_model="",
|
|
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="",
|
|
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="log_damp", min_error=3e-3)
|
|
|
|
|
|
# Helper functions to save and restore simulator state
|
|
def save_state(device):
|
|
state = {}
|
|
for region in ["Silicon", "Oxide", "Molding"]:
|
|
state[region] = {
|
|
"Potential": list(devsim.get_node_model_values(device=device, region=region, name="Potential"))
|
|
}
|
|
state["Silicon"]["Electrons"] = list(devsim.get_node_model_values(device=device, region="Silicon", name="Electrons"))
|
|
state["Silicon"]["Holes"] = list(devsim.get_node_model_values(device=device, region="Silicon", name="Holes"))
|
|
return state
|
|
|
|
def restore_state(device, state):
|
|
for region in ["Silicon", "Oxide", "Molding"]:
|
|
devsim.set_node_values(device=device, region=region, name="Potential", values=state[region]["Potential"])
|
|
devsim.set_node_values(device=device, region="Silicon", name="Electrons", values=state["Silicon"]["Electrons"])
|
|
devsim.set_node_values(device=device, region="Silicon", name="Holes", values=state["Silicon"]["Holes"])
|
|
|
|
def robust_solve(max_iters=40):
|
|
# Stage 1: Pre-conditioning solve (loose tolerance)
|
|
try:
|
|
devsim.solve(type="dc", absolute_error=1e13, relative_error=1e-1, charge_error=1e13, maximum_iterations=10)
|
|
except devsim.error:
|
|
pass
|
|
# Stage 2: Precision Newton solve
|
|
devsim.solve(type="dc", absolute_error=1e12, relative_error=3e-2, charge_error=1e13, maximum_iterations=max_iters)
|
|
|
|
# Solve zero-bias Drift-Diffusion state
|
|
print(f"[{temp_C} C] Solving initial zero-bias Drift-Diffusion...")
|
|
robust_solve(max_iters=40)
|
|
print(f"[{temp_C} C] Initial Drift-Diffusion converged.")
|
|
|
|
zero_bias_state = save_state(device)
|
|
|
|
# --- 5. Phase 1: Ramp Vds to target voltage at Vgs = 0.0V ---
|
|
vds_target = args.vds
|
|
print(f"[{temp_C} C] Ramping Vds from 0.0V to {vds_target:.3f}V (Vgs = 0.0V)...")
|
|
vds_curr = 0.0
|
|
vds_step = 0.1
|
|
|
|
# Sync gate contacts at 0V
|
|
for c in ["gate_Ox", "gate_Mold"]:
|
|
devsim.set_parameter(device=device, name=f"{c}_bias", value=0.0)
|
|
|
|
state_vds = save_state(device)
|
|
|
|
while vds_curr < vds_target:
|
|
vds_next = min(vds_curr + vds_step, vds_target)
|
|
for c in ["contR_Si", "contR_Mold"]:
|
|
devsim.set_parameter(device=device, name=f"{c}_bias", value=vds_next)
|
|
|
|
try:
|
|
robust_solve(max_iters=40)
|
|
vds_curr = vds_next
|
|
state_vds = save_state(device)
|
|
vds_step = min(vds_step * 1.5, 0.2)
|
|
except devsim.error:
|
|
# Revert and shrink step
|
|
for c in ["contR_Si", "contR_Mold"]:
|
|
devsim.set_parameter(device=device, name=f"{c}_bias", value=vds_curr)
|
|
restore_state(device, state_vds)
|
|
vds_step *= 0.5
|
|
if vds_step < 1e-4:
|
|
print(f"[{temp_C} C] Error: Vds ramping failed to converge at {vds_next:.3f}V. Aborting.")
|
|
sys.exit(1)
|
|
|
|
print(f"[{temp_C} C] Vds successfully ramped to 1.0V.")
|
|
|
|
# --- 6. Phase 2: Sweep Vgs from 0.0V to 10.0V at constant Vds = 1.0V ---
|
|
print(f"[{temp_C} C] Sweeping Vgs from 0.0V to 10.0V at constant Vds = 1.0V...")
|
|
vgs_target = 10.0
|
|
vgs_curr = 0.0
|
|
vgs_step = 0.05 # Start with fine step to capture the subthreshold region precisely
|
|
|
|
vgs_list = [vgs_curr]
|
|
# Extract current at Vgs = 0
|
|
i_n = devsim.get_contact_current(device=device, contact="contR_Si", equation="ElectronContinuityEquation")
|
|
i_p = devsim.get_contact_current(device=device, contact="contR_Si", equation="HoleContinuityEquation")
|
|
ids_list = [i_n + i_p]
|
|
|
|
g_av = devsim.get_edge_model_values(device=device, region="Silicon", name="AvalancheGeneration")
|
|
edge_vol = devsim.get_edge_model_values(device=device, region="Silicon", name="EdgeNodeVolume")
|
|
iav_list = [2.0 * sum(g * v for g, v in zip(g_av, edge_vol))]
|
|
|
|
state_vgs = save_state(device)
|
|
|
|
while vgs_curr < vgs_target:
|
|
vgs_next = min(vgs_curr + vgs_step, vgs_target)
|
|
for c in ["gate_Ox", "gate_Mold"]:
|
|
devsim.set_parameter(device=device, name=f"{c}_bias", value=vgs_next)
|
|
|
|
try:
|
|
robust_solve(max_iters=40)
|
|
vgs_curr = vgs_next
|
|
state_vgs = save_state(device)
|
|
|
|
# Extract current
|
|
i_n = devsim.get_contact_current(device=device, contact="contR_Si", equation="ElectronContinuityEquation")
|
|
i_p = devsim.get_contact_current(device=device, contact="contR_Si", equation="HoleContinuityEquation")
|
|
total_curr = i_n + i_p
|
|
|
|
# Integrated avalanche current
|
|
g_av = devsim.get_edge_model_values(device=device, region="Silicon", name="AvalancheGeneration")
|
|
edge_vol = devsim.get_edge_model_values(device=device, region="Silicon", name="EdgeNodeVolume")
|
|
total_av_curr = 2.0 * sum(g * v for g, v in zip(g_av, edge_vol))
|
|
|
|
vgs_list.append(vgs_curr)
|
|
ids_list.append(total_curr)
|
|
iav_list.append(total_av_curr)
|
|
|
|
print(f" Vgs = {vgs_curr:.3f}V, Ids = {total_curr:.4e}A, I_av = {total_av_curr:.4e}A (step: {vgs_step:.3f}V)")
|
|
|
|
# Adaptive step size
|
|
vgs_step = min(vgs_step * 1.5, 0.2)
|
|
gc.collect()
|
|
|
|
except devsim.error:
|
|
# Revert and shrink step
|
|
for c in ["gate_Ox", "gate_Mold"]:
|
|
devsim.set_parameter(device=device, name=f"{c}_bias", value=vgs_curr)
|
|
restore_state(device, state_vgs)
|
|
vgs_step *= 0.3
|
|
if vgs_step < 1e-4:
|
|
print(f"[{temp_C} C] Error: Vgs sweep failed to converge at Vgs = {vgs_next:.3f}V. Aborting.")
|
|
break
|
|
|
|
# Save CSV output for this temperature
|
|
csv_file = os.path.join(OUT_DIR, f"transfer_vgs_sweep_{temp_C:+.1f}C.csv")
|
|
np.savetxt(csv_file, np.column_stack((vgs_list, ids_list, iav_list)),
|
|
header="Vgs(V),Ids(A),Iav_eval(A)", delimiter=",", comments="")
|
|
print(f"[{temp_C} C] Saved transfer data to {csv_file}")
|
|
|
|
# Save final field visualization for ParaView
|
|
visual_variables = {"x", "y", "Potential", "Electrons", "Holes", "NetDoping", "Emag", "logEmag", "Jmag", "logJmag"}
|
|
for reg in ["Silicon"]:
|
|
devsim.element_from_edge_model(edge_model=opts['Jn'], device=device, region=reg)
|
|
devsim.element_from_edge_model(edge_model=opts['Jp'], device=device, region=reg)
|
|
devsim.element_model(device=device, region=reg, name="Jn_x", equation=f"{opts['Jn']}_x")
|
|
devsim.element_model(device=device, region=reg, name="Jn_y", equation=f"{opts['Jn']}_y")
|
|
devsim.element_model(device=device, region=reg, name="Jp_x", equation=f"{opts['Jp']}_x")
|
|
devsim.element_model(device=device, region=reg, name="Jp_y", equation=f"{opts['Jp']}_y")
|
|
devsim.element_model(device=device, region=reg, name="Jmag", equation="((Jn_x + Jp_x)^2 + (Jn_y + Jp_y)^2)^(0.5)")
|
|
devsim.element_model(device=device, region=reg, name="logJmag", equation="log(Jmag + 1e-20) / log(10.0)")
|
|
|
|
for reg in ["Oxide", "Molding"]:
|
|
devsim.element_model(device=device, region=reg, name="Jn_x", equation="0.0")
|
|
devsim.element_model(device=device, region=reg, name="Jn_y", equation="0.0")
|
|
devsim.element_model(device=device, region=reg, name="Jp_x", equation="0.0")
|
|
devsim.element_model(device=device, region=reg, name="Jp_y", equation="0.0")
|
|
devsim.element_model(device=device, region=reg, name="Jmag", equation="0.0")
|
|
devsim.element_model(device=device, region=reg, name="logJmag", equation="-20.0")
|
|
|
|
for reg in ["Silicon", "Oxide", "Molding"]:
|
|
devsim.element_from_edge_model(edge_model="EField", device=device, region=reg)
|
|
devsim.element_model(device=device, region=reg, name="Emag", equation="(EField_x^2 + EField_y^2)^(0.5)")
|
|
devsim.element_model(device=device, region=reg, name="logEmag", equation="log(Emag + 1e-20) / log(10.0)")
|
|
|
|
tec_file = os.path.join(OUT_DIR, f"transfer_sweep_{temp_C:+.1f}C_final.tec")
|
|
devsim.write_devices(file=tec_file, type="tecplot", include_test=lambda x: x in visual_variables)
|
|
print(f"[{temp_C} C] Saved 2D field visualization to {tec_file}")
|
|
|
|
print(f"[{temp_C} C] Transfer sweep finished successfully!")
|