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