import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
from scipy.integrate import trapz
from scipy.interpolate import interp1d
import os
# =============================================================
# 物理常数
# =============================================================
# ε0: 真空介电常数 (F/cm)
epsilon_0 = 8.854e-14
# q: 电子电荷 (C)
q = 1.602e-19
# k: 玻尔兹曼常数 (J/K)
k = 1.380649e-23
# T: 温度 (K)
T = 300
# β = q/(kT): 热电压的倒数 (V^{-1})
beta = q / (k * T)
# =============================================================
# 材料参数 (根据文件2.pdf更新)
# =============================================================
params = {
'N_D': 1e16, # cm^{-3}, 衬底施主掺杂浓度 (公式4.4)
't_F': 200e-7, # cm, 铁电层厚度 (200 nm)
'epsilon_F': 200, # 铁电层相对介电常数
'epsilon_i': 3.9, # 绝缘层相对介电常数 (公式4.4)
'epsilon_Si': 11.7, # Si相对介电常数
#'P_s': 10e-6, # C/cm^2, 饱和极化
#'P_r': 8.5e-6, # C/cm^2, 剩余极化
#'E_c': 50e3, # V/cm, 矫顽电场
#'delta': 75e3, # V/cm, 铁电材料参数
'EOT': 4e-7, # cm, 等效氧化层厚度
'A_f': 100e-8, # cm^2, 栅面积
'W': 10e-4, # cm, 沟道宽度
'L': 1e-4, # cm, 沟道长度
'mu': 200, # cm^2/V·s, 载流子迁移率
'n_i': 1.5e10, # cm^{-3}, 本征载流子浓度
}
# =============================================================
# 计算派生参数
# =============================================================
# 绝缘层电容 C_i = A_f * ε_i * ε_0 / EOT (公式4.7)
params['C_i'] = params['A_f'] * epsilon_0 * params['epsilon_i'] / params['EOT']
# 德拜长度 L_D = sqrt(ε_Si ε_0 φ_t / (q N_D)) (公式4.4)
phi_t = k * T / q # 热电压
params['L_D'] = np.sqrt(epsilon_0 * params['epsilon_Si'] * phi_t / (q * params['N_D']))
# 费米势 ψ_b = -φ_t * ln(N_D / n_i) (公式4.10)
params['psi_b'] = -phi_t * np.log(params['N_D'] / params['n_i'])
# =============================================================
# 读取外部P-E关系数据
# =============================================================
pe_file = r"D:\code\phi\pe.txt"
if not os.path.exists(pe_file):
raise FileNotFoundError(f"P-E关系文件未找到: {pe_file}")
# 加载数据
pe_data = np.loadtxt(pe_file)
E_ext = pe_data[:, 0] # MV/m
P_ext = pe_data[:, 1] # C/m²
# 将数据拆分为上扫和下扫分支
max_idx = np.argmax(E_ext)
E_up = E_ext[:max_idx+1] # 上扫分支
P_up = P_ext[:max_idx+1]
E_down = E_ext[max_idx:] # 下扫分支
P_down = P_ext[max_idx:]
# 转换为代码使用的单位
E_up = E_up *28.05* 1e4 # MV/m -> V/cm
P_up = P_up * 0.26* 1e-4 # C/m² -> C/cm²
E_down = E_down *28.05*1e4
P_down = P_down * 0.26 *1e-4
# 创建插值函数
pe_interp_up = interp1d(E_up, P_up, kind='linear', fill_value="extrapolate")
pe_interp_down = interp1d(E_down, P_down, kind='linear', fill_value="extrapolate")
# =============================================================
# 外部极化函数
# =============================================================
def external_polarization(E, direction):
if direction == 'up':
return pe_interp_up(E)
else: # down
return pe_interp_down(E)
# =============================================================
# C-V 特性仿真 (对应文件2.pdf第4.1节)
# =============================================================
def Q_s(psi_s):
"""
半导体表面电荷密度 (C/cm^2)
根据文件公式(4.4):
Q_s(ψ_s) = -sign(ψ_s) * [√2 * ε0 * ε_i] / (β * L_D) *
sqrt[ (n_i²/N_D²)(e^{-βψ_s} + βψ_s - 1) + (e^{βψ_s} - βψ_s - 1) ]
"""
term = (params['n_i']**2 / params['N_D']**2) * (np.exp(-beta * psi_s) + beta * psi_s - 1) + \
(np.exp(beta * psi_s) - beta * psi_s - 1)
return -np.sign(psi_s) * (np.sqrt(2) * epsilon_0 * params['epsilon_i']) / (beta * params['L_D']) * np.sqrt(term)
def cv_simulation(V_g, polarization_func, direction='up'):
"""
C-V特性仿真
基于文件公式(4.1)-(4.8)
"""
C_total = np.zeros_like(V_g)
psi_s0, E_f0 = 0.1, 0.1 # 初始猜测值
for i, V_g_val in enumerate(V_g):
def equations(vars):
psi_s, E_f = vars
# 公式(4.3): P(E_f) = -Q_s(ψ_s)
P = polarization_func(E_f, direction)
Q_s_val = Q_s(psi_s)
eq1 = P + Q_s_val
# 公式(4.2): V_F = E_f * t_F
V_f = E_f * params['t_F']
# 公式(4.3): V_I = -Q_s(ψ_s) * EOT / (ε_0 ε_i)
V_i = -Q_s_val * params['EOT'] / (epsilon_0 * params['epsilon_i'])
# 公式(4.1): V_G = ψ_s + V_I + V_F
eq2 = V_g_val - (psi_s + V_i + V_f)
return [eq1, eq2]
# =============================================================
# 改进的牛顿法求解器
# =============================================================
def newton_solver(equations, initial_guess, max_iter=50, tol=1e-6):
"""
牛顿法求解非线性方程组
equations: 函数,返回方程残差 [f1, f2]
initial_guess: 初始猜测 [x0, y0]
返回: 解 [x, y]
"""
x, y = initial_guess
for i in range(max_iter):
# 计算函数值
f = np.array(equations([x, y]))
# 检查收敛
if np.linalg.norm(f) < tol:
return x, y
# 数值微分计算雅可比矩阵
dx = 1e-6
J = np.zeros((2, 2))
# df1/dx
f1_x = equations([x + dx, y])[0]
f1 = equations([x, y])[0]
J[0, 0] = (f1_x - f1) / dx
# df1/dy
f1_y = equations([x, y + dx])[0]
J[0, 1] = (f1_y - f1) / dx
# df2/dx
f2_x = equations([x + dx, y])[1]
f2 = equations([x, y])[1]
J[1, 0] = (f2_x - f2) / dx
# df2/dy
f2_y = equations([x, y + dx])[1]
J[1, 1] = (f2_y - f2) / dx
# 求解线性方程组 J * delta = -f
try:
delta = np.linalg.solve(J, -f)
except np.linalg.LinAlgError:
# 奇异矩阵,使用伪逆
delta = np.linalg.pinv(J) @ (-f)
# 更新解
x += delta[0]
y += delta[1]
# 应用物理约束
x = max(min(x, 20.0), -20.0) # 表面势约束
y = max(min(y, 1e7), -1e7) # 铁电电场约束
# print(i,np.linalg.norm(f))
# 未收敛
raise RuntimeError(f"牛顿法未在 {max_iter} 次迭代内收敛")
# =============================================================
# 改进的 C-V 特性仿真
# =============================================================
def cv_simulation(V_g, polarization_func, direction='up'):
C_total = np.zeros_like(V_g)
# 存储前一个解用于初始猜测
prev_psi_s, prev_E_f = 0.1, 0.1
for i, V_g_val in enumerate(V_g):
def equations(vars):
psi_s, E_f = vars
P = polarization_func(E_f, direction)
Q_s_val = Q_s(psi_s)
V_f = E_f * params['t_F']
V_i = -Q_s_val * params['EOT'] / (epsilon_0 * params['epsilon_i'])
eq1 = P + Q_s_val
eq2 = V_g_val - (psi_s + V_i + V_f)
return [eq1, eq2]
try:
# 基于栅压的智能初始猜测
if i == 0:
# 第一点:根据栅压符号猜测
if V_g_val > 0:
initial_guess = [0.5, 1] # 正栅压时表面势和电场较大
else:
initial_guess = [-0.5, 1] # 负栅压时表面势和电场较小
else:
# 使用前一个解作为初始猜测
initial_guess = [prev_psi_s, prev_E_f]
# 使用牛顿法求解
psi_s, E_f = newton_solver(equations, initial_guess)
prev_psi_s, prev_E_f = psi_s, E_f
# 使用公式(4.5)计算半导体电容
term = (params['n_i']**2 / params['N_D']**2) * (-np.exp(-beta * psi_s) + 1) + \
(np.exp(beta * psi_s) - 1)
# 添加数值保护
if term < 0:
term = 1e-6 # 避免负值
C_si = (params['A_f'] * epsilon_0 * params['epsilon_i']) / (np.sqrt(2) * params['L_D']) * np.sqrt(term)
# 铁电层电容
C_f = params['A_f'] * epsilon_0 * params['epsilon_F'] / params['t_F']
# 绝缘层电容
C_i = params['C_i']
# 计算总电容
C_total[i] = 1 / (1/C_f + 1/C_i + 1/C_si)
except Exception as e:
if i > 0:
C_total[i] = C_total[i-1]
else:
C_total[i] = 1e-12
print(f"CV求解失败: V_g={V_g_val}, 错误: {e}")
return C_total
# =============================================================
# 改进的 I-V 特性仿真
# =============================================================
def solve_surface_potential(V_g, V, polarization_func, direction):
def equations(vars):
psi_s, E_f = vars
P = polarization_func(E_f, direction)
Q_s_val = Q_s(psi_s)
V_i = -Q_s_val * params['EOT'] / (epsilon_0 * params['epsilon_i'])
V_f = E_f * params['t_F']
eq1 = P + Q_s_val
eq2 = V_g - (psi_s + V_i + V_f)
return [eq1, eq2]
# 基于栅压的智能初始猜测
if V_g > 0:
initial_guess = [0.5, 1] # 正栅压
else:
initial_guess = [-0.5, -1] # 负栅压
try:
psi_s, E_f = newton_solver(equations, initial_guess)
return psi_s
except:
# 回退到物理合理的默认值
return abs(params['psi_b'])
def pao_sah_integrand(psi, V):
"""
改进的Pao-Sah被积函数,添加数值保护
"""
# 计算公式(21)中的内部表达式
term1 = np.exp(beta * psi) - beta * psi - 1
term2 = (params['n_i']**2 / params['N_D']**2) * np.exp(beta * V) * \
(np.exp(-beta * psi) + beta * psi * np.exp(-beta * V) - 1)
# 添加数值保护
total_term = term1 + term2
# print(total_term)
if total_term <= 0:
# 物理上应为正,处理数值误差
total_term = 1e-10
# 计算电场 ξ(ψ,V)
prefactor = (2 * params['N_D'] * k * T) / (epsilon_0 * params['epsilon_Si'])
xi = np.sqrt(prefactor * total_term)
# 公式(20)被积函数
numerator = (params['n_i']**2 / params['N_D']) * np.exp(-beta * (psi - V))
# 避免除零
if xi < 1e-10:
return 0
return numerator / xi
def pao_sah_integral(V_ds, V_g, polarization_func, direction):
"""
改进的Pao-Sah双积分,使用自适应积分
"""
# 使用自适应步长
num_V_points = max(10, int(V_ds * 10)) # 步长随V_ds自适应
V_points = np.linspace(0, V_ds, num_V_points)
integral_val = 0
for i in range(len(V_points)):
V_val = V_points[i]
try:
psi_s = solve_surface_potential(V_g, V_val, polarization_func, direction)
except:
psi_s = abs(params['psi_b'])
if psi_s > params['psi_b']:
# 自适应步长:在变化大的区域增加点数
num_psi_points = max(20, int((psi_s - params['psi_b']) * 100))
psi_points = np.linspace(params['psi_b'], psi_s, num_psi_points)
#print(psi_s)
# 计算被积函数值
integrand_vals = []
for psi in psi_points:
val = pao_sah_integrand(psi, V_val)
# 处理异常值
if not np.isfinite(val) or val < 0:
val = 0
integrand_vals.append(val)
inner_integral = trapz(integrand_vals, psi_points)
# 计算dV (梯形法)
if i == 0:
dV = V_points[1] - V_points[0]
elif i == len(V_points) - 1:
dV = V_points[-1] - V_points[-2]
else:
dV = (V_points[i+1] - V_points[i-1]) / 2
integral_val += inner_integral * dV
return integral_val
def iv_simulation(V_g, V_ds, polarization_func, direction='up'):
"""
I-V特性仿真
基于文件公式(4.9):
I_D = q μ (W/L) * [Pao-Sah双积分]
"""
I_d = np.zeros_like(V_ds)
for i, V_ds_val in enumerate(V_ds):
try:
integral = pao_sah_integral(V_ds_val, V_g, polarization_func, direction)
I_d[i] = q * params['mu'] * (params['W'] / params['L']) * integral
except Exception as e:
if i > 0:
I_d[i] = I_d[i-1]
else:
I_d[i] = 1e-12
print(f"IV求解失败: V_ds={V_ds_val}, 错误: {e}")
return I_d
# =============================================================
# 主程序
# =============================================================
if __name__ == "__main__":
# 设置扫描电压范围
V_g_cv = np.linspace(-1.6, 0.4, 30) # C-V栅压扫描
V_ds_iv = np.linspace(-0.5, 0, 20) # I-V漏压扫描
V_g_iv = 3.0 # I-V固定栅压
V_ds_transfer = 0.1 # 转移特性固定漏压
print("开始C-V仿真 (上扫)...")
C_up = cv_simulation(V_g_cv, external_polarization, direction='up')
print("开始C-V仿真 (下扫)...")
C_down = cv_simulation(V_g_cv[::-1], external_polarization, direction='down')
C_down = C_down[::-1] # 反转结果以匹配电压顺序
print("开始I-V仿真 (ON状态)...")
I_d_on = iv_simulation(V_g_iv, V_ds_iv, external_polarization, direction='up')
print("开始I-V仿真 (OFF状态)...")
I_d_off = iv_simulation(V_g_iv, V_ds_iv, external_polarization, direction='down')
print("开始转移特性仿真...")
V_g_transfer = np.linspace(0, 5, 20)
I_d_transfer = np.zeros_like(V_g_transfer)
for i, V_g_val in enumerate(V_g_transfer):
I_d_transfer[i] = iv_simulation(V_g_val, [V_ds_transfer], external_polarization, direction='up')[0]
# 绘制结果
plt.figure(figsize=(15, 12))
# 1. P-E曲线(新增子图)
plt.subplot(2, 3, 1)
plt.plot(E_up, P_up, 'b-', label='Up Sweep') # 直接使用原始数据
plt.plot(E_down, P_down, 'r--', label='Down Sweep')
plt.xlabel('Electric Field (MV/m)') # 单位已正确
plt.ylabel('Polarization (C/m²)') # 单位已正确
plt.title('Ferroelectric P-E Curve')
plt.grid(True)
plt.legend()
# 2. C-V曲线(原1号子图调整到2号位置)
plt.subplot(2, 3, 2)
plt.plot(V_g_cv, C_up * 1e6, 'b-',label='Up Sweep')
plt.plot(V_g_cv, C_down * 1e6, 'r--', label='Down Sweep')
plt.xlabel('Gate Voltage (V)')
plt.ylabel('Capacitance (μF)')
plt.title('MFIS Capacitor C-V Characteristics')
plt.grid(True)
plt.legend()
# 3. I-V曲线 (ON/OFF状态)(原2号子图调整到3号位置)
plt.subplot(2, 3, 3)
plt.plot(V_ds_iv, I_d_on * 1e6, 'g-', label='ON State')
plt.plot(V_ds_iv, I_d_off * 1e6, 'k--', label='OFF State')
plt.xlabel('Drain Voltage (V)')
plt.ylabel('Drain Current (μA)')
plt.title(f'FeMFET I-V Characteristics (V_g = {V_g_iv}V)')
plt.grid(True)
plt.legend()
# 4. 不同栅压下的I-V曲线(原3号子图调整到4号位置)
plt.subplot(2, 3, 4)
V_g_list = [1.0, 2.0, 3.0, 4.0]
for V_g_val in V_g_list:
I_d = iv_simulation(V_g_val, V_ds_iv, external_polarization, direction='up')
plt.plot(V_ds_iv, I_d * 1e6, label=f'V_g = {V_g_val}V')
plt.xlabel('Drain Voltage (V)')
plt.ylabel('Drain Current (μA)')
plt.title('FeMFET I-V at Different Gate Voltages')
plt.grid(True)
plt.legend()
# 5. 转移特性曲线(原4号子图调整到5号位置)
#plt.subplot(2, 3, 5)
#plt.semilogy(V_g_transfer, np.abs(I_d_transfer) * 1e6, 'mo-')
#plt.xlabel('Gate Voltage (V)')
#plt.ylabel('Drain Current (μA)')
#plt.title(f'Transfer Characteristics (V_ds = {V_ds_transfer}V)')
#plt.grid(True)
plt.tight_layout()
# 保存图像
output_dir = r"D:\code\phi"
if not os.path.exists(output_dir):
os.makedirs(output_dir)
plt.savefig(os.path.join(output_dir, "mfisfet_simulation_results.png"), dpi=300)
print(f"结果已保存至: {os.path.join(output_dir, 'mfisfet_simulation_results.png')}")
plt.show()
# 保存数据
np.savetxt(os.path.join(output_dir, "cv_data.txt"),
np.column_stack((V_g_cv, C_up * 1e6, C_down * 1e6)),
header="V_g(V) C_up(μF) C_down(μF)")
np.savetxt(os.path.join(output_dir, "iv_data.txt"),
np.column_stack((V_ds_iv, I_d_on * 1e6, I_d_off * 1e6)),
header="V_ds(V) I_on(μA) I_off(μA)")
print("仿真数据已保存")现在求解器有巨大问题,积分和解方程都不收敛。能不能用matlab写一个。matlab有没有能直接解方程的工具