1、原理：载流子浓度=载流子总数的可能最大值载流子出现概率，参考这里
2、计算细节，先对材料算结构优化，然后做单点能计算输出电荷密度，最后做非自洽计算输出DOS
非自洽计算输出DOS看这里
Global Parameters
ISTART = 1 (Read existing wavefunction, if there)
ISPIN = 1 (Non-Spin polarised DFT)
ICHARG = 11 (Non-self-consistent: GGA/LDA band structures)
LREAL = Auto (Projection operators: automatic)
ENCUT = 450 (Cut-off energy for plane wave basis set, in eV)
PREC = Accurate (Precision level: Normal or Accurate, set Accurate when perform structure lattice relaxation calculation)
LWAVE = .F. (Write WAVECAR or not)
LCHARG = .F. (Write CHGCAR or not)
ADDGRID= .TRUE. (Increase grid, helps GGA convergence)
# LVTOT = .TRUE. (Write total electrostatic potential into LOCPOT or not)
# LVHAR = .TRUE. (Write ionic + Hartree electrostatic potential into LOCPOT or not)
# NELECT = (No. of electrons: charged cells, be careful)
# LPLANE = .TRUE. (Real space distribution, supercells)
# NWRITE = 2 (Medium-level output)
# KPAR = 2 (Divides k-grid into separate groups)
# NGXF = 300 (FFT grid mesh density for nice charge/potential plots)
# NGYF = 300 (FFT grid mesh density for nice charge/potential plots)
# NGZF = 300 (FFT grid mesh density for nice charge/potential plots)
NCORE = 4
IVDW = 11
Electronic Relaxation
ISMEAR = -5 (Gaussian smearing, metals:1)
SIGMA = 0.05 (Smearing value in eV, metals:0.2)
NELM = 200 (Max electronic SCF steps)
NELMIN = 6 (Min electronic SCF steps)
EDIFF = 1E-06 (SCF energy convergence, in eV)
# GGA = PS (PBEsol exchange-correlation)
Ionic Relaxation
NSW = 0 (Max ionic steps)
IBRION = -1 (Algorithm: 0-MD, 1-Quasi-New, 2-CG)
ISIF = 3 (Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions)
EDIFFG = -0.02 (Ionic convergence, eV/AA)
# ISYM = 2 (Symmetry: 0=none, 2=GGA, 3=hybrids)
NEDOS = 1500
#LORBIT = 11
#LSCALAPACK = .FALSE.
IOPTCELL = 1 1 0 1 1 0 0 0 0
#NELECT=92.8 #94
#AMIX = 0.05
#BMIX = 0.0001
2.1 载流子计算脚本
##载流子计算脚本20251204
import numpy as np
kb = 8.617333e-5 # eV K⁻¹
T = 300 # K
S = 17.3*4 # 单胞面积 Ų，改为自己数值
EF=-0.0881 #费米能级
EC=0-EF #导带底相对费米能级
EV=3.1000-EF #价带顶底相对费米能级
dos= np.loadtxt(‘DOSCAR’, skiprows=6)
E = dos[:,0] – EF # 能量轴，E=0 为 费米能级
tdos= dos[:,1] # 总 DOS
de = E[1]-E[0] # DOSCAR 能量步长
f = lambda e: 1/(np.exp(e/kb/T)+1)
n = np.trapz(tdos[EC>=0],f(E[EC>=0]))/S*1e16
p =np.trapz(tdos[EV<=0]*(1-f(E[EV<=0])))/S*1e16
print(‘n = {:.2e} cm⁻²’.format(n))
print(‘p = {:.2e} cm⁻²’.format(p))