1994 PISCES-2ET and Its Application Subsystem
1994 PISCES-2ET and Its Application Subsystem
1994 PISCES-2ET and Its Application Subsystem
Application Subsystems
Table of Contents
5 Numerical Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.1 High and Zero Frequency AC Analyses . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.1.1 General Principle of AC Small Signal Analysis . . . . . . . . . . . . . . . . . . . . . 43
5.1.2 Matrix Transformation to Improve Diagonal Dominance . . . . . . . . . . . . . 46
5.1.3 Pre-conditioned TFQMR Method for AC Analysis . . . . . . . . . . . . . . . . . . 47
5.1.4 Zero Frequency AC Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.2 Discretization Scheme for Carrier and Energy Fluxes . . . . . . . . . . . . . . . 50
5.3 Normalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.4 Newton Projection Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.5 Global Data Structure in PISCES-2ET . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6 Simulation Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6.1 Substrate Current of MOS Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6.2 SOI Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
6.3 Cylindrically Symmetric LED . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
6.4 AlInAs/GaInAs MODFET . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
DOPING . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .117
ELECTRODE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .127
ELIMINATE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .130
END (QUIT). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .132
EXTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .133
IMPACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .136
INCLUDE. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .138
INTERFACE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .139
LOAD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .141
LOG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .144
MATERIAL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .146
MESH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .151
METHOD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .157
MOBILITY. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .164
MODELS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .170
OPTIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .180
PLOT.1D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .183
PLOT.2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .192
PRINT. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .196
REGION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .200
REGRID . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .206
SOLVE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .212
SPREAD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .222
SYMBOLIC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .226
TITLE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .230
VECTOR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .231
X.MESH, Y.MESH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .233
Reference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263
16 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .291
16.1 Introduction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .291
16.2 Single Stage CMOS Inverter (DC Analysis) . . . . . . . . . . . . . . . . . . . . . .291
16.3 High Frequency GaAs MESFET (AC Analysis) . . . . . . . . . . . . . . . . . . .292
16.4 SRAM Cell (Transient Analysis) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .294
16.5 GaAs/AlGaAs LED (Transient Analysis). . . . . . . . . . . . . . . . . . . . . . . . .298
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .315
PISCES-2ET is a dual energy transport (for carrier temperatures and lattice thermal diffusion) version
of PISCES-II developed based on Stanford’s 9009 version and Intel’s enhanced 8830 version. This
work is completed under the support of SRC (Semiconductor Research Corp.), ARO (Army Research
Office), and ARPA (Advanced Research Projects Agency) and with close collaboration from Intel and
HP. There are many new features available in the 2ET version, predominately, the capabilities to
simulate the carrier and lattice temperatures and heterostructures in compound semiconductors. Hence,
various non-stationary phenomena such as hot carrier effects and velocity overshoot can be analyzed
using this program. The electrical behavior of optoelectronic devices can also be simulated with
reasonable accuracy. Most of the material parameters have been calibrated and thoroughly surveyed
with the help from industry. Other new features which can be found in PISCES-2ET include:
This document is organized as follows. CHAPTER 2 discusses the theoretical background for
DUET transport model, which is based on the moment approach to solving Boltzmann Transport
Equation (BTE) and is applicable to heterostructure analysis also. In CHAPTER 3, physical models
used in the program are described. The emphasis is on the various mobility models, many parameters
of which have been calibrated by industrial users. Material parameters for four common compound
semiconductors, namely GexSi1-x, AlxGa1-xAs, AlxIn1-xAs, and GaxIn1-xAsyP1-y, used in
heterostructure and optoelectronic devices are discussed in CHAPTER 4. Discussion of carrier
transport in heterostructures and issues related to Fermi-Dirac (FD) statistics are delayed until
APPENDIX A to provide a complete picture of DUET model yet not to inundate the concise forms in
the main text with complex coefficients required by using the FD distribution. In CHAPTER 5,
relevant numerical techniques are presented, including the algorithms for small signal AC analyses.
Besides, the details in parameterization of global arrays used in the device solver are provided. Finally
four simulation examples are included in CHAPTER 6. The first two examples compare the simulated
results to the experimental data regarding the carrier temperature effects in FET structures to establish
the confidence level of the energy transport model used in the code. The other two examples are about
the simulation of heterostructures including one application in a realistic, cylindrically symmetric LED
(light emitting diode) structure. Two short appendices, one dealing with the mathematical properties
of Fermi integral (APPENDIX B) and the other listing several expressions referred to in this document
from the previous PISCES reports (APPENDIX C), are included for users’ reference. Finally a user’s
manual provides the detailed syntax and usage of all commands available in the program.
Authors are grateful to many contributions from both academia and industry during the
development of PISECE-2ET. The DUET transport model has been developed in collaboration with
Prof. K. Hess’ group in University of Illinois at Urbana. One of the key developers, Dr. E. Kan, who
recently joined our group at Stanford, made numerous contributions to the code including helping
parameterize the data structure and model calibration. Drs. T. Thurgate of Intel and M. Tan of HP have
been working closely with us from the very beginning of this project and provided access to their codes
and experimental data. Various industrial connections under the sponsorship of SRC are greatly
appreciated, including, but by no means being complete, Drs. R. Lowther of Harris Semiconductor, I.
Lim of Motorola, and K. L. Chen of OCD in HP. Dr. K. Wu of this group at Stanford developed the
algorithm for high frequency small signal AC analysis, and Mr. A. Mujtaba helped review the section
of mobility models. Their contributions are greatly appreciated. Finally, the continuous funding from
SRC through contract #SRC 93-SJ-116 is one of the key factors to the success of this program
development and the grant from ARO through contract DAAL 03-91-G-0152 plays a critical role in
extending the simulation capabilities from silicon to compound materials and heterostructure. The
GaAs work was supported through ARPA contract DAAL 01-91-K-0145. Their support are gratefully
acknowledged.
The DUET model, a carrier transport model in semiconductors, is developed based on the moment
approach to solving Boltzmann Transport Equation (BTE). It uses six state variables to describe the
status of a semiconductor device. These six variables are: electrostatic potential, ψ, carrier
concentrations, n and p, carrier temperatures, Tn and Tp, and lattice temperature, TL, and they are
functions of space and time. All other device characteristics such as terminal I-V characteristics and
circuit model parameters can be calculated from the knowledge of the distribution of these basic
variables. To determine the distribution of these variables under applied bias, six independent
equations are required together with proper boundary conditions. It is well established that with the
drift-diffusion (DD) carrier transport model, Shockley semiconductor equations, i.e., Poisson’s
equation and carrier continuity equations, govern the distribution of ψ, n, and p. The carrier
concentrations can also be replaced, equivalently, by their respective quasi-Fermi levels, φ n and φ p ,
in classical distribution (either Boltzmann or Fermi-Dirac) functions. With the temperatures for both
carriers and lattice introduced as independent variables, three more equations are needed and they can
be derived from the energy balance principle. In this chapter, we will first introduce two (kinetic)
energy balance equations for carriers and one thermal diffusion equation for the lattice. Various macro
quantities such as carrier and energy densities will be defined using the distribution function at the
quasi-thermal equilibrium. The auxiliary expressions for transport are provided. Finally the proper
boundary conditions for solving differential equations are discussed.
h ∂f0
f ( r, k ) = f 0 ( r, k ) – τ ( r, k ) ------------*-k ⋅ ∇ f 0 – ------------*- -------- k ⋅ E
h
(2.1)
2πm 2πm ∂ε
where r and k are spatial (position) and wavenumber vectors, respectively, f 0 is the distribution
function at the quasi-thermal equilibrium in both the real ( r ) and momentum ( k ) space, τ is the
relaxation time which depends on both r and k , and so is defined in a microscopic sense, E is the
electric field and ε is the carrier kinetic energy, e.g. for electrons, ε = E – E C where E C is the
conduction band edge. Note that in this manual, we use bold faced E ( E ) and it components to
represent the electric field, while the plain faced E is used for the total energy of carriers. All other
symbols in the above expression have conventional meanings.
The quasi-thermal equilibrium function, f 0, is chosen from one of the classical distribution
functions and for the present we will use the Boltzmann distribution function. In APPENDIX A, we
will describe how to extend the distribution function to Fermi-Dirac statistics. The Boltzmann
distribution function has the following form in r and k space:
E – E Fn ε – ( E Fn – E C )
f 0 ( r, k ) = 2 exp – ------------------ = 2 exp – ----------------------------------- (2.2)
kBT n kBT n
where E Fn and T n are the quasi-Fermi energy and temperature for electrons, respectively. In our later
analysis, it is often desirable to express f 0 in terms of the carrier concentration. We now proceed this
transformation. From definition, it is easy to conclude that the carrier concentration is the zero-
moment of the distribution function in k space. Thus taking electrons as an example, we have
E Fn – E C
n = N C exp --------------------- (2.3)
kBT n
where N C is the effective density of states for the conduction band and is found to be
3⁄2
2πm n ( T L )k B T n
*
NC = 2 ------------------------------------
- (2.4)
2
h
where T L is the lattice temperature and because in our DUET model, T L is not necessarily the same
*
as T n , we have explicitly indicated the lattice temperature dependence of the effective mass, m n .
ε
f 0 ( r, k ) = f 0 ( n ( r ), T n ( r ), T L ( r ), ε ) = ---------------------------- exp – ------------
2n
(2.5)
N C ( T n, T L ) k B T n
This form of distribution function constitutes the basis for our model discussion.
∂w n
--------- = – ∇ ⋅ s n + j n ⋅ E n – u wn (2.6)
∂t
where w is the kinetic energy density and s is the energy flux, and both can be defined and evaluated
from the carrier distribution function. The last two terms on the right hand side (RHS) of the above
equation represent the energy conversion and net loss rates, respectively. The familiar Joule heat term,
j n ⋅ E n where E n is the electric field acting on electrons1, actually represents the rate of conversion
from the electrostatic potential to the kinetic energy, and u w is the rate of net loss (loss minus
1. In heterostructures, the electric field acting on carriers is in general different for electrons and
holes.
generation) for kinetic energy due to carrier recombination/generation and energy exchange with the
lattice through phonon scattering. The first step in the development of carrier transport model is then
to define and determine the expressions for various quantities used in the above balance equation.
DUET model is based on the Stratton’s description of distribution function at non-thermal equilibrium
(Eq. (2.1) and refer to [1] for details). The concept of the carrier temperature can directly be deduced
from the definition of the kinetic energy density in this approach, and it turns out that the parameter
Tc, where c stands for either n or p, used in the distribution function, f 0,at quasi-thermal equilibrium
(Eqs. (2.2) and (2.5)) is the carrier temperature following the ideal gas kinetics. The carrier
concentration (n), kinetic energy density (wn), current density (jn), and energy flux (sn) can all be
evaluated from the distribution function, f in Eq. (2.1), from their respective definitions by integration
in the momentum (i.e. k ) space:
E Fn – E C
∫fd
3
n = p = N C exp --------------------- (2.7)
kBT n
2
p 3
∫
3
wn = ---------*- fd p = --- nk B T n (2.8)
2m n 2
q
∫
3
j n = – -----*- pf d p
mn
3 * T (2.9)
= qD n ∇n – --- qnD n ∇lnm n + nµ n ∇E C + qnD n ∇T n
2
= nµ n ∇E Fn + qnµ n Q n ∇T n
2
1 p
∫
3
s n = -----*- p ---------*- f d p = – P n T n j n – κ n ∇T n (2.10)
m n 2m n
where we have introduced the carrier momentum, p = ( h ⁄ 2π )k , and the integral element volume,
3 3
d p = ( dk x dk y dk z ) ⁄ ( 2π ) and constant effective mass in (crystal) momentum space is assumed1.
T
The transport coefficients D (diffusion constant), µ (carrier mobility), D (thermal diffusivity or called
“Soret” coefficient in [1]), Q (thermopower), P (thermoelectric power), and κ (thermal conductivity)
1. This assumption is not required in the DUET model (see [2]) but is used for present implementa-
tion of the code.
are all defined and can be evaluated from the even part of the distribution function f in the momentum
space ( f 0 in Eq. (2.1) for DUET model), and are thus related to each other. For continuity of
presentation, the expressions for these coefficients are not given here, rather they can be found in
APPENDIX A. We now consider the rates for energy generation g w and loss r w , from which
uw = r w – gw (2.11)
r w is related to both carrier recombination and energy transferring from carriers to the lattice. At
present, three recombination mechanisms are considered in the model and they are Shockley-Reed-
Hall (SRH), Auger, and radiative recombinations. For g w , only the impact ionization, which is the
inverse process of Auger recombination, is taken into account1. With all relevant mechanisms
identified, we are able to list the complete set of equations for the DUET model and their auxiliary
expressions as follows:
Poisson’s equation:
+ –
∇ ⋅ ( – ε∇ψ ) = q ( p – n + N D – N A ) (2.12)
where dielectric constant, ε , is not to be confused with carrier kinetic energy, ε .
∂n 1
----- = --- ∇ ⋅ j n – u (2.13)
∂t q
∂p 1
------ = – --- ∇ ⋅ j p – u (2.14)
∂t q
where u is the net recombination rate of electron-hole pairs.
∂w n
--------- = – ∇ ⋅ s n + j n ⋅ E n – u wn (2.15)
∂t
1. Thus, photo-generation as discussed in Section 3.3.2 applies only to the conservation of carrier
population.
∂w p
--------- = – ∇ ⋅ s p + j p ⋅ E p – u wp (2.16)
∂t
∂T L 3 3
c L --------- = ∇ ⋅ ( κ L ∇T L ) + u SRH --- k B T n + E g ( T L ) + --- k B T p +
∂t 2 2
(2.17)
wn ( T n ) – wn ( T L ) w p ( T p ) – w p ( T L )
----------------------------------------- + ------------------------------------------
τ wn τ wp
where τ wn and τ wp are energy relaxation times for electrons and holes, respectively.
The following auxiliary expressions are needed for current and energy flux densities, and
energy loss rate in Eqs. (2.15)-(2.17).
Current densities:
3 * T
j n = qD n ∇n – --- qnD n ∇lnm n + nµ n ∇E C + qnD n ∇T n
2 (2.18)
= nµ n ∇E Fn + qnµ n Q n ∇T n
3 * T
j p = – q D p ∇p + --- qpD p ∇lnm p + pµ p ∇E V – qpD p ∇T p
2 (2.19)
= pµ p ∇E Fp – qpµ p Q p ∇T p
s n = – P n T n j n – κ n ∇T n (2.20)
s p = P p T p j p – κ p ∇T p (2.21)
where D and µ are related by Einstein relationship (see Eq. (A.23) for a general expression) and Q is
also linked to µ in a more complex way (Eq. (A.43)), other coefficients have more simple expressions
if only Boltzmann statistics is considered.
3 kB
P n = P p = P = --- ----- (2.22)
2q
κ n = nk B T n µ n P (2.23)
κ p = pk B T p µ p P (2.24)
It is apparent that the carrier mobility, µ, is the key parameter to other transport coefficients, and the
mobility can be obtained either from the theoretical calculation or more often through empirical or
semi-empirical formulation by fitting the analysis/simulation results to the experimental data. We will
discuss in a great detail the mobility modeling in CHAPTER 3. Finally, the rate of net loss of carrier
kinetic energy is modeled as follows:
3 3
u wn = --- ( u SRH + u rad )k B T n – ( u n, Auger – g n, imp ) E g ( T L ) + --- k B T p –
2 2
(2.25)
3 wn ( T n ) – wn ( T L )
--- g k T – ----------------------------------------
-
2 p, imp B n τ wn
3 3
u wp = --- ( u SRH + u rad )k B T p – ( u p, Auger – g p, imp ) E g ( T L ) + --- k B T n –
2 2
(2.26)
3 w p(T p) – w p(T L)
--- g k B T p – -----------------------------------------
-
2 n, imp τ wp
where u SRH , u rad , and u Auger are net carrier loss rates due to SRH, Auger, and radiative
recombinations, respectively, g imp is the generation rate due to the impact ionization (II). Note that for
both Auger recombination and II, we need to distinguish if they are caused by electrons or holes. The
last term on RHS of each above expression represents the energy exchange between the carriers and
lattice.
A special case for the above general description (Eqs. (2.15)-(2.17)) is when only the lattice
temperature is of concern, i.e., the carriers can be considered to have the same temperature as the
lattice. Such a scenario is often encountered in, say, the simulation of power devices. We can then lump
Eqs. (2.15)-(2.17) by requiring all temperatures being the same (designated T) and need to solve the
following thermal diffusion equation only.
∂
----- c L + --- ( n + p )k B T = ∇ ⋅ ( κ L ∇T – s n – s p ) + j n ⋅ E n + j p ⋅ E p + u SRH E g ( T )
3
(2.27)
∂t 2
Together with Poisson’s (Eq. (2.12)) and carrier continuity equations (Eqs. (2.13)-(2.14)), these four
equations are solved for four state variables: ψ, n, p, and T (for both lattice and carriers).
Normally the Dirichlet thermal boundary condition, i.e., the known carrier and lattice
temperatures, is applied to both thermal and electrical contacts. In many cases, however, the heat
conduction through the intimate contacts to the device is modeled by a thermal resistance. In this case,
the lattice (and carrier) temperature at the contact becomes unknown. The relationship between the
heat flux and temperatures at the contact (T L, cont ) and thermo-reservoir (the environment temperature,
T env ) is specified as follows:
T env – T L, cont
F th = --------------------------------
- (2.28)
R th
where the direction of the heat flux, F th , is defined as going into the device. In the present model, the
carrier temperature is always assumed the same as the lattice temperature at the contact. There may be
two types of contacts to the device: electrical and thermal ones. An electrical contact must be a thermal
contact also, but the contrary is not necessarily true. For those boundaries of a device region where
there is no electrical/thermal contact specified, a reflective boundary condition is assumed as is the
case for electrical boundary conditions.
µ ( N , T L, E ⊥, E || /T c ) = f ( µ 0 ( N , T L, E ⊥ ), E || /T c ) (3.1)
where E || and E ⊥ are the longitudinal and transverse components of the electric field with respect to
(w.r.t) the current direction, µ0, which depends on N, TL, and E ⊥ , is called the low field mobility
because when E || → 0 , µ → µ 0 , Tc is the carrier temperature with the subscript c representing either
n or p for electrons and holes, respectively, and the symbol E || ⁄ T c indicates the dependency is either
on E || or on Tc but not on both. In some models such as ones developed by Intel, the above expression
is simplified to the following form to separate mobility reductions due to different field components:
In the following, we will first consider the low field mobility models without taking into account the
effect of E ⊥ . We then consider the effect of the transverse field dependence and finally discuss the
longitudinal field or carrier-temperature dependence.
TABLE 3.1
µn (cm2 V -1 s-1) µp (cm2 V -1 s-1)
Silicon / GexSi1-x 1500 450
GaAs / AlxGa1-x As 8500 400
InAs / AlxIn1-xAs 22600 250
InP / GaxIn1-xAsyP1-y 4500 150
be used instead with this option. The range of the data for doping concentrations in silicon is from
14 21
N = 1 ×10 to 1 ×10 cm –3 . If the doping concentration falls out of the above range then the closest
boundary value for the mobility is used, otherwise the following interpolation scheme is applied.
µ0 ( N 2 ) – µ0 ( N 1 ) N
µ 0 ( N ) = µ 0 ( N 1 ) + ----------------------------------------- log ------ (3.3)
log ( N 2 ⁄ N 1 ) N1
where N 2 > N > N 1. This model is valid for silicon only. Note that if any of the following parameters:
ccsmob, analytic, arora, and user1, is used with or without conmob specified in the model
card, that model takes precedence and the precedency of these models is user1, arora, analytic,
and ccsmob with each preceding one overriding the trailing ones.
y – y int – d
γ = 0.5 ⋅ erfc -------------------------- (3.5)
λ
where the silicon and SiO2 interface is assumed to be parallel to the x-axis, y int is the y coordinate of
the interface, d and λ are two parameters: one for offset and the other for the characteristic length.
Their default values are d = 0.06µ and λ = 0.03µ , and both can be accessed by users via parameters
int.off and char.int in mobility card, respectively.
The surface and bulk mobilities have the same form of dependence on N and TL as follows:
d
TL α( N )
µ ( N , T L ) = ----------------------µ max + ---------------------- µ min (3.6)
1 + α( N ) 1 + α( N )
c b
where α ( N ) = ( N ⁄ N 0 ) T L and hereinafter T L = T L (in K) ⁄ 300 , the normalized lattice
temperature. But N has different meaning for bulk and surface. For bulk, N is simply the total doping
concentration, while for surface N must include the effect of the interface charge in the following
manner:
N = N + cQ int (3.7)
-2) as specified by parameter qf in the interface card,
where Q int is interface charge (in units of cm
and c is a scaling factor in units of cm-1 with default value of 107 and can be changed through
parameter qss.conc in mobility card. Note that the effect of the surface mobility is taken into
account only for those nodes which are directly (in the y-direction) under the Si/SiO2 interface. There
is one exception, however. That is, when the flag intelmob.par is set in model card, the surface
mobility (µ srf in Eq. (3.4)) uses a different expression from Eq. (3.6). The new expression has the form
of
–a
µ srf = µ 0 T L (3.8)
where parameters µ 0 and a for silicon are listed in Table 3.2.
Table 3.2
µ 0 (cm2 V -1 s-1) a
electrons 783.5 1.67
holes 247.4 1.13
There are six parameters in Eq. (3.6): µ max , µ min , N 0 , b, c, and d, and all are dependent upon
the carrier type (electrons or holes) and region (bulk or surface). The parameter values used in the code
for silicon are listed in Table 3.3.
Table 3.3
µ dlt
µ 0 ( N , T L ) = µ min + --------------------------------
-
α
(3.9)
1 + (N ⁄ N 0)
b
where parameters µ min , µ dlt , N0, and α are all functions of TL in the form of aT L where both a and b
are constants. This model is later extended to apply to GaAs by Yu [6] based on the available measured
data at 77 and 300 K. The parameter values for silicon and GaAs are listed in Table 3.4.
Table 3.4
† †
µ min µ dlt N0 (cm-3) α
– 0.57 – 2.33 17 2.4 – 0.146
electrons 88T L 1252T L 1.25 ×10 TL 0.88T L
– 0.57 – 2.23 17 2.4 – 0.146
silicon holes 54.3T L 407T L 2.35 ×10 T L 0.88T L
– 0.7475 – 2.687 16 3.535 – 0.1441
electrons 2136T L 6331T L 7.345 ×10 T L 0.6273T L
– 1.124 – 2.366 17 3.690
GaAs holes 21.48T L 331.2T L 5.136 ×10 T L 0.8057
† units of cm2 V -1 s-1
1.025
µ ( N , T L, n, p ) = µ L ( T L ) -----------------------------------------------------------------------------
1.43
- – 0.025 (3.10)
1 + { ( X ( N , T L, n, p ) ) ⁄ 1.68 }
–α
µ L ( T L ) = µ L0 T L (3.11)
6µ L [ µ I ( N , T L ) + µ ccs ( T L, n, p ) ]
X ( N , T L, n, p ) = ------------------------------------------------------------------------------
- (3.12)
µ I ( N , T L )µ ccs ( T L, n, p )
3⁄2 2 2 –1
AT L BT BT L
µ I ( N , T L ) = --------------- ln 1 + ---------L- – --------------------
- (3.13)
N N N + BT
2
L
17 3⁄2
2 ×10 T L 8 2 –1 ⁄ 3 –1
µ ccs ( T L, n, p ) = ---------------------------- [ ln { 1 + 8.28 ×10 T L ( pn ) }] (3.14)
pn
where the parameters µ L0 , α , A, and B are listed in Table 3.5.
Table 3.5
Note that in the program, N is taken as the min (N, 1019) where N is in units of cm-3.
two such models available in the program, both were developed by Intel (intelmob and
intelmob.par in model card). We then consider the Lombardi model (lombardi) which shares
the same feature as Intel’s ones in that both transverse and longitudinal field effects are considered
simultaneously in the code. In all these three models, the transverse field can be computed either based
on the structure or by definition. We will elaborate this at the end of this section.
–β
r perp = ( 1 + E ⊥ ⁄ E crit ) (3.15)
4 4
where Ecrit is a parameter with value of 4.2×10 V/cm for electrons and 3×10 V/cm for holes in
silicon, and β = 0.5 . These parameters are all accessible to users through mobility card.
The second model specified by intelmob.par has different expression for µ srf (Eq. (3.8))
and more user-accessible parameters for the reduction due to the transverse field:
1
r perp = ---------------------------------------α- (3.16)
1 + ( E ⊥ ⁄ E univ )
where the parameters α and Euniv are listed in Table 3.6.
Table 3.6
α E univ ( V/cm )
5
electrons 1.02 5.71 ×10
5
holes 0.95 2.57 ×10
Both Intel’s mobility models use µ 0 ( N , T L, E ⊥ ) = r perp ( E ⊥ )µ 0 ( N , T L ) for the low-field mobility.
1 1 1 1
-------------------------------- = ------------------------------------ + --------------------- + ------------------------ (3.17)
µ ( N , T L, E ⊥ ) µ ac ( N , T L, E ⊥ ) µ srf ( E ⊥ ) µ 0 ( N , T L )
where µac is the mobility due to the acoustic phonon scattering, which depends on both the transverse
field and lattice temperature, µsrf is the one due to the surface scattering and is the function of the
transverse field only, and µ0(N, TL) can use any model in Section 3.1.1. The expressions for µac and
µsrf are as follows
β
B αN
µ ac ( N , T L, E ⊥ ) = ------ + ----------------- (3.18)
E⊥ T E1 ⁄ 3
L ⊥
δ
µ srf ( E ⊥ ) = ------2 (3.19)
E⊥
The parameters used in the above expressions are listed in Table 3.7 for electrons and holes,
respectively.
Table 3.7
B α β δ
7 5 14
electrons 4.75 ×10 1.74 ×10 0.125 5.82 ×10
7 5 14
holes 9.93 ×10 8.84 ×10 0.0317 2.05 ×10
Note that the units for N, TL and E ⊥ in the above two equations and in conjunction with the use of
parameter values in Table 3.7 are cm-3, K, and V/cm, respectively.
There are now problems remain unspecified in the above Intel’s and Lombardi models. That
is the evaluation of the transverse field and the reduction due to the longitudinal field. We discuss these
issues in the following subsections.
E× j
E ⊥ = ---------------- (3.20)
j
PISCES-2ET uses this formula to evaluate E ⊥ in the above three mobility models if the parameter
strfld is not specified in model card. Such evaluation is computationally costly and might also
cause numerically instability. Another alternative and simplification is to evaluate the normal field
component based on the geometric proximity to the Si/SiO2 interface. Thus when strfld is
specified, the following formula is used:
– ( y – y int ) ⁄ L
E⊥ = E ye (3.21)
where y is the coordinate in the axis perpendicular to the interface, E y is the transverse field at the
interface (yint), and L is the characteristic length which has a value of 0.5µ.
1
E ⊥, eff = ----- ( ζQ d + ηQ i ) (3.22)
ε si
where Qd is the charge density (per unit surface area) in the inversion layer due to the ionized
impurities and Qi is the charge density due to the inverted (minority) carriers, and both ζ and η are
physics-based fitting parameters. It is usually taken that ζ = 1 for both carriers and η = 0.5 for
electrons and η = 0.33 for holes. The dependence of the mobility on E ⊥, eff is modeled in the
following way by using Mathiessen’s rule:
1 1 1 1
------------- = -------- + ------- + --------
µ 0, inv µ ph µ sr µ im
6 α1 6 α2
1 1 ×10 1 1 ×10
= ------------ -------------- + ------------ -------------- (3.23)
µ ref 1 E ⊥, eff µ ref 2 E ⊥, eff
18 – 1 12 α3
1 1 ×10 1 ×10
+ ------------ ---------------- ----------------
µ ref 3 N B N i
where µph is the mobility due to the phonon scattering, µsr is the one due to the surface roughness
scattering, and µim is due to the scattering of the charged impurities in the inversion layer. And NB is
the doping density in the substrate (units of cm-3), and Ni is the inverted carrier density per unit surface
area in the channel (units of cm-2). The other parameters are listed in Table 3.8.
Table 3.8
The modeling of mobility dependence on the transverse field can generally be categorized as
two classes: local or non-local field dependence. The local field dependence refers to the fact that the
transverse electric field is computed locally, whereas non-local field dependence is to find an effective
transverse field to be used in the mobility formulation and this effective field is computed based on the
knowledge over the entire inversion (or accumulation) layer in the silicon substrate at the Si/SiO2
interface. All the models discussed so far are of local field nature. In next section we will discuss two
non-local models, both developed at University of Texas, Austin.
µ 0 ( N , T L, E ⊥ ) E⊥ – E0 dµ
µ ( N , T L, E ⊥, E || ) = --------------------------------- - ---------0-
- + --------------------------------------- (3.24)
µ 0 E || µ E 2 3 ⁄ 2 dE
1 + -----------
2 ⊥
1 + -----------
0 ||
v sat v sat
where the dependence of µ 0 is listed only once. Note that all four parameters, N , T L, E ⊥ , and E || are
functions of space, so are µ and µ 0 . The non-locality is represented by the quantity of E 0 which is
defined as the transverse electric field at the edge (i.e. bottom) of the inversion layer. The inversion
layer edge is in turn determined by the criterion that the inverted carrier concentration starts to drop
below the doping density. The expression for µ 0 is
µ 0 ( N , T L, E ⊥ ) = ------------------------
- + 3.2 ×10 --- T L
1 –9 p 1 ⁄ 2 –1
(3.25)
–5 ⁄ 2 z
1150T L
where
3⁄2 –8 N f
p = 0.09T L + 1.5 ×10 ----------------
1⁄4
- (3.26)
n TL
–5
0.039T 1.24 ×10
z = -------------------L- + ------------------------------
- (3.27)
E ⊥ + E 0 E ⊥ + E 01 ⁄ 3
------------------- -------------------
2 2
where n is the inverted carrier density in the inversion layer and in this particular case is the electron
concentration for the n-channel MOSFET, and Nf is the fixed interface charge per unit area at the Si/
SiO2 interface. To use this model, users need to provide the location of the gate oxide through
parameters ox.left, ox.right, and ox.bottom in the model card.
–1
µ 0 ( N , T L, E ⊥ ) = -------- + ------- + ------
1 1 1
(3.28)
µ ph µ sr µ C
where µph is the mobility due to the scattering by phonons in both bulk and surface and fixed interface
charge, and µsr is the one due to the surface roughness scattering, and µC is to the screened Coulomb
scattering. The detailed expressions for those three mobilities can be found in [13].
µ ( N , T L, E ⊥, E || ) = r para ( E || )µ 0 ( N , T L, E ⊥ ) (3.29)
and the reduction factor r para is modeled as (refer to [14])
µ 0 ( N , T L, E ⊥ )E || β –1 ⁄ β
r para = 1 + ---------------------------------------- (3.30)
v sat
where in silicon β is 1.395 for electrons and 1.215 for holes, and the saturation velocity, v sat , has the
expression in Table 3.9. The parameter β can be changed by the user through two parameters,
b.electrons and b.holes, in the model card.
v sat E || 4
µ 0 ( N , T L ) + -------- ------
E E 0
µ ( N , T L, E || ) = ----------------------------------------------------- (3.31)
E || 4
1 + ------
E 0
7 4
where E0 is the critical field and has a default value of 4kV/cm, and v sat = 1.13 ×10 – 1.2 ×10 T L . It
can be shown that when E > E 0 Eq. (3.31) leads to a negative differential mobility (NDM). This is the
default field dependent mobility model for electrons in GaAs and the value of E0 can be altered by the
user through the parameter E0 in the model card.
One problem related to the above formulation is that when applied to the simulation of GaAs
MESFET, the drain output characteristics (current vs. voltage) may exhibit an unrealistic zig-zag
behavior. One possible reason is that the correct mobility model in the bulk may not be suitable to the
carriers in the surface channel. In the program, we provide another mobility model for electrons in
GaAs, in which the carrier velocity approaches the saturation velocity with increased field
monotonically in a manner of the hyperbolic tangent function. The model formulation is as follows
[16] and does not exhibit NDM behavior.
v sat µ 0 ( N , T L )E ||
µ ( N , T L, E || ) = -------- tanh ------------------------------ (3.32)
E || v sat
This model can be invoked through parameter hypertang in the model card.
µ 0 ( N , T L, E ⊥ )
µ ( N , T L, E ⊥, T c ) = --------------------------------------------------------------------------
- (3.33)
3
1 + --- α ( N , T L, E ⊥ )k B ( T c – T L )
2
where the subscript c stands for carriers, α is defined as being proportional to µ0 in the following way:
µ 0 ( N , T L, E ⊥ )
α ( N , T L, E ⊥ ) = ---------------------------------
2
- (3.34)
qτ w v sat ( T L )
where τw is the energy relaxation time and vsat is the saturation velocity for carriers, respectively. Their
default values for silicon are listed in Table 3.9 where T L is in units K.
Table 3.9
τ w ( ps ) v sat ( cm/s )
6 3
electrons 0.42 8.64 ×10 – 2.68 ×10 T L
6 3
holes 0.25 10.65 ×10 – 2.68 ×10 T L
2
µ ( N , T L, E ⊥, E || ) = µ 0 ( N , T L, E ⊥ ) -------------------------------------------------------------------------- (3.35)
µ 0 ( N , T L, E ⊥ )E || 2
1 + 1 + 4 ----------------------------------------
v sat
–α
κL( T L ) = κ0T L (3.36)
where κ0 is the thermal conductivity at TL= 300 K and varies from material to material while α is a
constant (1.2) for all materials. Table 3.10 lists κL for different materials as implemented in the code.
Table 3.10
G = α n nv n + α p pv p
(3.37)
where vn and vp are the electron and hole velocities, respectively, and the impact ionization rates, αn
and αp, which are defined as the number of electron-hole pairs generated by a carrier per unit distance
traveled. The II rate thus defined is a strong function of the local electric field or more accurately (as
most researchers now agree) of local carrier temperature. PISCES-2ET provides several II rate models
having either local field or carrier temperature dependence. In the following we will first discuss the
more classical local field dependent model and provide the relevant material parameters and then the
carrier temperature dependent models.
Table 3.11
1. Correction has been made to the original formulation (Eq. (80) on p. 59 of [18]).
2. We do not consider in Eq. (3.38) the situation of j ⋅ E < 0 , i.e. the current is against the field direc-
tion.
Table 3.11
5 kB T c – T L
E eff = --- ----- ------------------ (3.39)
2 q γ Lw
where L w = τ w v sat is called the energy relaxation length and γ is an adjustable parameter (for fitting
the experimental data) which can be accessed by the user using impjt.ratio in the model card
(default: 1.0). The second way is to approximate the actual carrier (drift) velocity with their saturation
velocity and to use carrier-temperature dependent ionization rate as follows:
–T c ⁄ Bc
α ( T c ) = Ac e (3.41)
where c stands for either n or p. The values of A and B are listed inTable 3.12. This model can be
Table 3.12
A (cm-1) B (K)
6 6
electrons 3.8 ×10 1.75 ×10
7 6
holes 2.25 ×10 3.26 ×10
The third modeling approach in developing carrier temperature dependent II is to abandon the
generation rate dependence on the carrier (average) velocity at all. The rationing behind this approach
is that when the carrier energy is higher than some threshold value the generation rate should be
determined by the number of high energy carriers only no matter which direction they travel. This
model can be invoked by parameter imp.nt in the model card, and the formulation and its
parameters are as follows:
B –C n ⁄ T n B –C p ⁄ T p
G = n An T n n e + p A pT p pe (3.42)
where the parameters A, B, and C are listed Table 3.13
Table 3.13
A B C
11
electrons 5.5 ×10 4.0 0.8
13
holes 2.0 ×10 4.0 1.0
3.3.2 Photo-Generation
This model is invoked by parameter photogen in the model card. The generation rate due to
photons is modeled as a forced generation term which is proportional to the photon flux, F ph , and the
intensity decays exponentially deep into the device along the incident path. The decay length is
characterized by the absorption coefficient, α. The current implementation assumes only y-
dependence:
– αy
g ( y ) = αF ph e (3.43)
-2 -1 -1
where Fph has units of cm s and α units of cm . Both Fph and α have to be specified by the user
through parameters flux and abs.coef in the model card.
where Et is the energy level for the recombination centers and Ei is the intrinsic Fermi Energy. Carrier
lifetimes, τn and τp, can be either constant or dependent upon the total doping concentration (parameter
consrh in model card). The formula for the doping dependent is as follows:
τ0
τ = -----------------------------------
- (3.45)
1 + N tot ⁄ N SRH
where Ntot is the total doping concentration, and the default values for parameters τ0 and Nsrh are the
same for both types of carriers in silicon and τ 0 = 1 ×10 s and N srh = 5 ×10 cm-3. If the interfacial
–7 16
recombination is to be included, it is through the change of carrier lifetime in the following way:
1 ds int 1
-------- = ----------
- + --- (3.46)
τ eff A τ
where A and d are partial area and interface length associated with the node in the element,
respectively, and s int is the interface recombination velocity.
2
u Aug ( n, p ) = ( c n n + c p p ) ( np – n i ) (3.47)
cm6 s-1,
– 31 – 32
where the default values for c n and c p in silicon are 2.8 ×10 and 9.9 ×10 respectively.
The minimum set of material parameters which have to be known in order to proceed the
device simulation include:
1. Electron affinity, χ , and conduction band edge offset due to the composition change.
2. Energy bandgap, E g .
* *
3. Effective masses for electrons and holes, m n and m p , from which the effective densities of states
for the conduction and valence bands, NC and NV, can be derived (Eq. (2.4)).
4. Static and high-frequency dielectric constants, ε 0 and ε ∞ .
5. Electron and hole mobilities, µ n and µ p , and their dependence on composition, doping density,
electric field, and temperature.
6. Minority carrier lifetime and corresponding coefficients for various recombination mechanisms
including Auger and radiative.
7. Coefficients for the impact ionization.
8. Saturation velocity, v sat , which is used as a parameter in certain field-dependent mobility models.
Before we address the detailed values and dependence of these parameters, however, we will
first discuss the interpolation scheme for determining the compound material parameters based on the
knowledge of the base materials which constitute the compound.
1
Q ( x, y ) = ---------------------------------------------- ⋅ {
x(1 – x) + y(1 – y)
x ( 1 – x ) [ ( 1 – y )T Ga x In1 – x P + yT Ga x In1 – x As ] + (4.1)
y ( 1 – y ) ( 1 – x ) [ ( 1 – x )T InAs y P1 – y + xT GaAs y P1 – y ] }
where GaxIn1-xAsyP1-y is used as an example for quaternary material. We now provide a number of
material parameters in PISCES-2ET as discussed below.
Table 4.1
2
αT L
E g ( T L ) = E g0 – ---------------
- (4.2)
β + TL
where TL is in K. When the data for the light and heavy holes are available, the overall hole effective
mass used in the program is determined by
*3 ⁄ 2 *3 ⁄ 2 *3 ⁄ 2
mp = m lh + m hh (4.3)
We have also defined three auxiliary materials: Ge, AlAs, and GaP, which are used to constitute the
desired ternaries and quaternary. Their material parameters relevant to simulation are listed below:
Table 4.2
Now we discuss how the parameters for the compound materials are determined, starting with
the bandgap.
1.16 – 0.945
E g ( x ) = E g, Si – ------------------------------ x x ≤ 0.245
0.245
0.945 – 0.87
= 0.945 – ------------------------------ ( x – 0.245 ) 0.245 < x ≤ 0.35
0.35 – 0.245
0.87 – 0.78
= 0.87 – --------------------------- ( x – 0.35 ) 0.35 < x ≤ 0.5
0.5 – 0.35
0.78 – 0.72 (4.4)
= 0.78 – --------------------------- ( x – 0.5 ) 0.5 < x ≤ 0.6
0.6 – 0.5
0.72 – 0.69
= 0.72 – --------------------------- ( x – 0.6 ) 0.6 < x ≤ 0.675
0.675 – 0.6
0.69 – 0.66
= 0.69 – --------------------------------- ( x – 0.675 ) 0.675 < x ≤ 0.735
0.735 – 0.675
= 0.66 x > 0.735
For bandgap of AlxIn1-xAs, as x is increased from 0 to 1 the band structure also undergoes a
transition from the direct to indirect bandgap and the critical composition is x = 0.56 . The expression
for Eg at room temperature is as follows [31]:
where when x < 0.56 the (direct) bandgap is determined by the Γ valley of the conduction band
( E g ( AlAs ) = 2.95eV , E g ( InAs ) = 0.36eV ), and when x > 0.56 the (indirect) bandgap is
determined by the X valley ( E g ( AlAs ) = 2.16eV , E g ( InAs ) = 1.37eV .)
Before giving the expression for the bandgap in the quaternary GaxIn1-xAsyP1-y. we first list
the expressions for the bandgaps of four ternary materials, on which GaxIn1-xAsyP1-y is based.
2
E g ( Ga x In 1 – x P ) = 1.35 + 0.643x + 0.786x (4.7)
2
E g ( Ga x In 1 – x As ) = 0.36 + 0.505x + 0.555x (4.8)
2
E g ( InAs y P 1 – y ) = 1.35 – 1.083 y + 0.091y (4.9)
2
E g ( GaAs y P 1 – y ) = 2.74 – 1.473 y + 0.146y (4.10)
And for the strained GaxIn1-xAs layer, the bandgap is evaluated according to the expression
provided by Corzine et al. [32]:
If either x or y in GaxIn1-xAsyP1-y has the value of zero or one, then one of the expressions,
Eqs. (4.7) - (4.10), is used to evaluate the bandgap, otherwise the following general formula is used to
evaluate the (lowest-direct) bandgap in GaxIn1-xAsyP1-y [26]:
1
E g ( x, y ) = ---------------------------------------------- { x ( 1 – x ) [ ( 1 – y )E g ( Ga x In 1 – x P ) + yE g ( Ga x In 1 – x As ) ]
x(1 – x) + y(1 – y) (4.12)
+ y ( 1 – y ) ( 1 – x ) [ ( 1 – x )E g ( InAs y P 1 – y ) + xE g ( GaAs y P 1 – y ) ] }
The values obtained from the above expression is for room temperature ( T L = 300K ).
usually assumed in this case that the amount of the edge shift for both the conduction and valence
bands are the same. For the current implementation, except of AlxGa1-xAs, AlxIn1-xAs, and GaxIn1-
xAsyP1-y, it is assumed that all bandgap change due to the composition change occurs only at the
valence edge due to the lack of sufficient experimental data. However, the bandgap narrowing due to
the heavy doping effect, is accounted for in the aforementioned manner for all materials, i.e., evenly
split. For AlxGa1-xAs / GaAs, ∆E C : ∆E V is taken as 0.6 : 0.4 and for AlxIn1-xAs / InAs, 0.706 : 0.294,
where ∆E is considered positive when the band edge is expanding, that is the bandgap becomes larger.
For GaxIn1-xAsyP1-y / InP the bandgap change is simply assumed evenly split among the conduction
and valence bands edges1. By designating the percentage of the shift of conduction band edge with
respect to the bandgap change as γ, we have the following relationship for the electron affinity as
function of the mole fractions:
χ ( x, y ) = χ ( 0, 0 ) – γ [ E g ( x, y ) – E g ( 0, 0 ) ] (4.13)
∆E C = – ∆χ (4.14)
∆E V = ∆χ + ∆E g (4.15)
*
m p( x)
--------------- = 0.48 + 0.31x (4.16)
m0
For AlxIn1-xAs / InAs system, we simply use the linear interpolation scheme for both
electrons and light and heavy holes. That is,
* * *
m ( x ) = ( 1 – x )m InAs + xm AlAs (4.17)
where the effective masses for InAs and AlAs can be found in Table 4.1 and Table 4.2, respectively.
For GaxIn1-xAsyP1-y / InP, the interpolation scheme based on the constituent binaries is used
for obtaining the effective masses for electrons and light and heavy holes. For electrons, we don’t
distinguish between the bulk and the strained layer. The interpolation scheme is as follows:
* * *
m n = ( 1 – x )ym n, InAs + ( 1 – x ) ( 1 – y )m n, InP +
(4.18)
* *
xym n, GaAs + x ( 1 – y )m n, GaP
For holes, the effective masses of the light and heavy holes in the bulk material are calculated
in the same manner. However, for the strained layer, they are taken as constants according to Corzine
in [35]1.
*
m lh = 0.078m 0 (4.19)
*
m hh = 0.165m 0 (4.20)
1. We have intentionally flipped the data for the light and hole holes because of suspicion of the mis-
take incurred in the cited paper.
ε1 – 1 ε2 – 1
1 + 2 x -------------- + ( 1 – x ) --------------
ε1 + 2 ε2 + 2
ε = ------------------------------------------------------------------------
- (4.21)
ε1 – 1 ε2 – 1
1 – x -------------- – ( 1 – x ) --------------
ε1 + 2 ε2 + 2
where ε1 is that of AlAs for both the materials, and ε2 is that of GaAs for AlGaAs and of InAs for
AlInAs. The formulation applies to both high and low frequencies.
The dielectric constants for GexSi1-x are obtained by using the linear interpolation for those
of Ge and Si. For the quaternary GaxIn1-xAsyP1-y, the interpolation scheme in Eq. (4.18) is used based
on the data for the constituent binaries.
4.6 Mobilities
Due to the lack of experimental data, the mobility dependences on composition, doping, lattice
temperature, and field are not complete for most of the compound materials. Whenever data are not
available constant mobility as listed in Section 3.1.1.1 is used. Readers are referred to [29] for data for
many other materials, even though the collection of data is rather out-of-date (1971).
In this chapter various numerical techniques and data structure will be discussed. Since the
development of PISCES has spanned over a decade, we only address here the improvement over the
previous versions (II and II-B). Interested readers are referred to reports [3] and [4].
F ( x, p ) = 0 (5.1)
where both the equation vector (also a column vector) F and variable vector x have dimension of N,
N
the number of all basic variables in the simulation region, i.e., F, x ∈ R , while the parameter vector
M
p has the dimension of M, no bigger than the total number of the device terminals, or p ∈ R . For
analyses involving the time evolution such as transient and AC analyses, the above equation must be
modified to include the time derivative part. Thus,
F+ f = 0 (5.2)
where
∂x
f = D ------ (5.3)
∂t
where D is a diagonal matrix (N× N) with rows corresponding to the Poisson’s equation zero and others
– 1 . We now apply bias of form
jωt
p = p + p̃e (5.4)
That is, an AC signal of angular frequency ω and amplitude p̃ (notice that the multi-dimensionality
represents the bias at all terminals) which is generally a vector of complex numbers, superimposed on
the DC bias, p . The response, that is the solution vector x, should also have the form of
jωt
x = x + x̃e (5.5)
where vector x is the DC response and are the real numbers, and x̃ is a phasor vector and are complex
numbers due to the change of phase in response to the input signal. For the small signal analysis,
p̃ → 0 where represents the norm of the vector, hence x̃ → 0 . Substituting Eqs. (5.4) and (5.5)
into Eq. (5.2) and applying the Taylor expansion to Eq. (5.2), by retaining only terms which have the
jωt
first order of e on obtains the following equation:
jωt
F ( x, p ) + [ F x ( x, p ) x̃ + F p ( x, p )∆p + jωD x̃ ]e = 0 (5.6)
Note that in the above we have replaced p̃ with a real number vector ∆p , i.e. we assume all the input
signals are in the same phase. Furthermore, because both x and p are vectors the derivatives of F x
and F p are actually matrices with F x being a square one while F p is generally not.
With the DC solution already known, i.e. F ( x, p ) = 0 , Eq. (5.6) can be simplified to the
following:
( F x + jωD ) x̃ = – F p ∆p (5.7)
and the only unknown variable is x̃ , the AC response. Because x̃ is a complex number, we can write
x̃ = x r + jx i where x r and x i are real and imaginary parts, respectively, and j = – 1 , Eq. (5.7) then
becomes a set of two equations in real number mathematics:
F x x r – ωD x i = – F p ∆p (5.8)
ωD x r + F x x i = 0 (5.9)
These two equations are coupled to each other so to find xr and xi, they have to be solved
simultaneously in principle. Because each of vectors xr and xi has dimension of N, then we must solve
a system of equations with dimension of 2N × 2N which needs to quadruple the memory storage as is
required for the DC analysis. Because the AC analysis is usually not performed as routinely as the DC
analysis, to set aside a large trunk of memory space solely for the AC analysis is not an economic
choice especially when the program is not written to take the advantage of dynamic memory allocation.
It is preferable, therefore, to use an iterative method to solve the above equation. One popular method,
also the one used in PISCES IIB is the so-called SOR (Successive Over-Relaxation) method [36]. The
essence of this method can better be seen from the following matrix equation form:
J – ωD x r B
Ax = = (5.10)
ωD J xi 0
where we have altered the symbols used in Eqs. (5.8)-(5.9) with the Jacobian matrix J = F x and the
stimulus B = – F p ∆p which has non-zero (real number) terms only for rows corresponding to the
Poisson’s equation. If the magnitude of the off-diagonal elements, i.e., ω D, in the above block matrix
equation are small compared to the diagonal elements, then the iterative approach by repeating the loop
of first solving
J x r = B + ωD x i (5.11)
for xr while fixing xi and then solving
J x i = – ωD x r (5.12)
for xi using the recently updated xr, will converge easily. But as the ω increases, the convergence
becomes poorer and eventually fails to converge when the frequency approaches the cutoff frequency,
fT, of the device under analysis [36]. A final note regarding SOR method is the use of relaxation factor.
Instead of employing the above simple iteration procedure (Eqs. (5.11) and (5.12)), the actual iteration
takes form of
k+1 k –1 k
xr = ( 1 – ν ) x r + νJ ( B + ωD x i ) (5.13)
k+1 k –1 k+1
xi = ( 1 – ν ) x i + νJ ( – ω D x r ) (5.14)
where ν is the relaxation factor and its value ranges from 0 to 2, and k is the iteration count. When
ν = 1 , Eqs. (5.13) and (5.14) are reduced to Eqs. (5.11) and (5.12). One of the most challenging tasks
in SOR method is to determine the optimized value of ν and we will discuss later other iteration
schemes which avoid using the relaxation factor at all.
I ωDD J
–1
T = –1 (5.15)
– ωDD J I
where I is the identity matrix and D is the diagonal matrix consisting of only the diagonal elements in
J
J. The resulting equation is thus
2 –1 –1
J + ω DD J D ωDD J J 0 x r B
AT x = = (5.16)
– ωDD –J 1 J 0 J + ω DD J D x i
2 –1
0
where J 0 = J – D J and the property of DB = 0 has been used to arrive at the above equation1. Note
–1
that because D J is a diagonal matrix so its inverse D J is trivial to compute. It becomes immediately
clear from Eq. (5.16) that although the magnitude of the off-diagonal blocks still increases
proportionally with ω, that of the diagonal entries may increase even faster due to the additional term
2
which is proportional to ω , thus avoiding the problem of diagonal blocks being overwhelmed by the
off-diagonals. There are three distinctive ranges of frequencies where the transformed matrix shows
different characteristics. At low frequencies, which can be quantified by all the non-zero entries in
–1
ωDD J being much less than unity, off-diagonal blocks can be considered as small perturbations with
respect to the diagonal blocks and iterative algorithms should do well with or without the
–1
transformation. At very high frequencies, where all the non-zero entries of ωDD J are much greater
than unity, the enhancement to the diagonal blocks makes the entire matrix become increasingly
diagonally dominant, which implies the improved performance for the iterative algorithm. At the
intermediate frequencies between these extremes, however, the situation is not as clearly cut since the
off-diagonal blocks can no longer be considered as small perturbations while the growth in the
enhancements to the diagonal blocks has not caught up with the growth in the off-diagonals.
1. The right hand side of Eq. (5.10) has been multiplied by T but the result remains the same as the
original.
–1
needed during iteration. Designating the pre-conditioner as P = M , then the matrix M for Eq. (5.16)
will be
J + ω 2 DD –J 1 D 0
M = 2 –1 (5.17)
0 J+ω DD J D
and
J 0
M = (5.18)
0 J
for Eq. (5.10). Readers can verify that by using these pre-conditioners, the resulting coefficient
matrices have identity matrices as their diagonal blocks. At the low frequency, the off-diagonal blocks
are very small in magnitude, leading to a quick convergence.
As an iterative algorithm, the SOR method has the advantage of simplicity and minimal
requirement for the data storage. In practice, however, we have found that a proper relaxation
coefficient is difficult to obtain. In contrast, there are several attractive advantages with the so-called
QMR (Quasi-Minimal-Residual) iterative algorithm [39]. Until quite recently, all existing iterative
algorithms either use non-minimization procedure in generating the update vectors, which results in a
very volatile residual reduction process, or use a full minimization procedure which requires
increasing data storage as the iteration proceeds. By applying a quasi-minimal residual approach,
QMR is able to retain the advantage of requiring only the fixed amount of data storage. At the same
time, the residual reduction process becomes very smooth with only occasional, insignificant
increases. Furthermore, like all other Krylov subspace methods based on the Lanzos vectors, no
external parameters are needed [39]. The particular QMR algorithm we have implemented in PISCES-
2ET is the TFQMR (Transpose-Free Quasi-Minimal-Residual) [38] with the pre-conditioner M in
either Eq. (5.17) or Eq. (5.18) applied before iteration. The detailed procedure of TFQMR with pre-
conditioner M is summarized as follows:
1. Initialization:
-1
(a) b = M B;
(b) y 1 = u 0 = r 0 = r̃ = b , x 0 = d 0 = 0 ;
(c) ρ 0 = τ 0 = r 0 , θ 0 = η 0 = 0 ;
2. For n = 1, 2, … do
–1
(a) ν n – 1 = M Ay 2n – 1 ;
(b) σ n – 1 = r̃ ⋅ v n – 1 , α n – 1 = ρ n – 1 ⁄ σ n – 1 ;
(c) y 2n = u n – 1 – α n – 1 v n – 1 ;
–1
(d) r n = r n – 1 – α n – 1 M A ( u n – 1 + y 2n – 1 ) ;
(e) ω 2n + 1 = r n , ω 2n = rn ⋅ rn – 1 ;
(f) For m = 2n – 1, 2n do
2
• d m = y m + ( θ m – 1 η m – 1 ⁄ α n – 1 )d m – 1 ;
2
• θ m = ω m + 1 ⁄ τ m – 1, c m = 1 ⁄ 1 + θ m ;
2
• τ m = ω m + 1 c m, η m = c m α n – 1 ;
• xm = xm – 1 + ηm d m ;
(g) ρ n = r̃ ⋅ r n, β n = ρ n ⁄ ρ n – 1 ;
(h) u n = r n + β n y 2n ;
(i) y 2n + 1 = u n + β n ( y 2n + β n y 2n – 1 ) ;
In [38], 1 + m τ m is given as an upper limit for the residual norm. In our implementation, the
residual norm is initially assumed to be equal to τ m (a coefficient of unity). When τ m drops below the
tolerance, the true residual norm is evaluated, which also gives out the real proportional coefficient
between the residual norm and τ m . If the residual norm is not below the tolerance, the iteration
continues with the new coefficient.
The effectiveness of the TFQMR with pre-conditioning in AC analysis has been tested in an
analysis of GaAs MESFET [40]. There is essentially no limitation as to how high the frequency can
be applied to the device during simulation and the agreement between the measured and simulated
results are very good.
F x C = –D xr (5.19)
Information such as the change of the storage charge from the quasi steady state analysis can
thus be obtained by using the zero-frequency AC analysis, which involves only two more solutions of
a linear equation with the coefficient matrix the same as the Jacobian for the DC solution but different
right hand side terms.
1
F n, 1 → 2 = ------- [ µ n2 B ( u n, 21 )n 2 T n2 – µ n1 B ( – u n, 21 )n 1 T n1 ] (5.20)
d 21
1
F p, 1 → 2 = ------- [ µ p1 B ( u p, 21 ) p 1 T p1 – µ p2 B ( – u p, 21 ) p 2 T p2 ] (5.21)
d 21
where F n = – j n ⁄ q and F p = j p ⁄ q, and B is the Bernoulli function defined as
x
B ( x ) = -------------
x
(5.22)
e –1
u c, 21 is the function of the potential difference and carrier temperature defined as follows:
ψ2 – ψ1
u c, 21 = ---------------------------------- (5.23)
( T c1 + T c2 ) ⁄ 2
where c stands for either n or p and d21 is the distance between nodes 1 and 2. The assumptions made
in the above discretization scheme are:
1. Mobility is constant along the grid line between nodes 1 and 2 (called as interval) and its value
is the average of the mobilities on two nodes.
2. E ( x ) ⁄ T c ( x ) ≈ const over the interval where x represents the position within the interval. In
particular,
E ψ1 – ψ2 2
----- = ------------------- ----------------------- (5.24)
Tc d 21 T c1 + T c2
3. Flux, F, is constant over the interval.
ψ1 + ψ2
s n, 1 → 2 = ------- µ n2 B ( u n, 21 )n 2 T n2 ------------------
- – c e T n2 –
1
d 21 2
(5.25)
ψ1 + ψ2
µ n1 B ( – u n, 21 )n 1 T n1 ------------------
- – c e T n1
2
ψ1 + ψ2
s p, 1 → 2 = ------- µ p1 B ( u p, 21 ) p 1 T p1 ------------------
- – c e T p1 –
1
d 21 2
(5.26)
ψ1 + ψ2
µ p2 B ( – u p, 21 ) p 2 T p2 ------------------
- – c e T p2
2
where c e = 5 ⁄ 2 – ν 1 and is set to a default value of 1.8 in the code. The form for the heat
flux:
κ
F heat, 1 → 2 = ------- ( T L1 – T L2 ) (5.27)
d 21
where the average κ over the interval is evaluated as follows:
–α
T L1 + T L2
κ L = κ 0 ------------------------ (5.28)
2 × 300
5.3 Normalization
This section briefly describes the normalization scheme used in PISCES-2ET. Internal variables
representing basic variables: ψ (as well as φ’s for quasi Fermi potentials), n, p, Tn, Tp, and TL, are all
normalized before they are used in the computation. Potentials are normalized by V th = k B T 0 ⁄ q
where T0 is the environment temperature which can be accessed through parameter temperature
in the model card. All temperatures are normalized using the environment temperature and carrier
concentrations are normalized using c scl = N C N V where NC and NV are for the base semiconductor
material2. The current density is normalized using the quantity of – qc scl V th . The reason for choosing
this normalization factor is related to the use of Scharfetter-Gummel scheme and minus sign is incurred
because the electron current density is used as the reference. Other quantities such as distance, time,
electric field, mobility are not normalized.
p = p + ∆p (5.29)
one can obtain the following approximate equation:
Because the solution at p has been found, i.e., F ( x, p ) = 0 , the increment of x due to p can
be found by solving the following system of equations:
F x ∆x = – F p ∆p (5.31)
This is a linear system of equations with the same coefficient matrix as for the DC solution
and F p is easy to find because F has few equations which are related to p. The above formulation can
also be applied to the low-frequency small signal analysis because when ∆ p → 0 , the ratio of ∆x ⁄ ∆p
can be found exactly by solving the above equation for ∆x. Foe more detailed discussion of this
method, readers are referred to [41].
-------------------------p2conf.h------------------------------
c There are four parameters that affect the memory usage of
c PISCES-2ET:
c AVGNB - average neighbors of each node point including
c itself. For a quad-based mesh with control level
The total memory usage during the run time of PISCES-2ET is approximately:
Note that since the full Newton iteration is the default method owing to its robustness, the
memory space required is proportional to the square of the number of equations and to the power 1.5
of the number of points in the worst case (when MAX1D is equal to the square root of the number of
points, i.e., the mesh points are evenly distributed in X and Y directions). The present setup (MAXPT
= 3000) will take about 24M-byte memory space. Since most of the memory allocation is shared in
a common block called “tmpco” (in various.h files), users can still run DUET model with the present
setup at least to about 1,200 points. Most of the modern workstations have more than 24M-byte main
memory, so the swapping during execution should be minimal. If the application cannot be bounded
by this setup, users need only to change the above four parameters and remake (recompile) the
PISCES-2ET executable. PISCES-2ET will automatically report any insufficient memory allocation
during the first encounter (usually at the symbolic factorization). Users can then modify the reported
dimensions and make further trade-off accordingly. If the application is always much smaller than the
present setup, we do not suggest to recompile. If MAXPT becomes too small and the factored LU
matrix in full Newton iteration is no longer dominant in the memory usage, users probably will have
some complaints from the compiler about array with negative dimensions. This is because the common
“tmpco” space is shared in many routines and the patch array TMPPAD is used to make every
common declaration to have the same length to avoid compiler errors. If users have to do this (probably
due to an old 8MB main memory machine) and the compiler indeed complains the negative dimension,
users need to change the last line in p2conf.h, where the size of the common block “tmpco” is
defined symbolically. Other constraints will have minimal impact on the memory usage. They are
listed below for reference.
---------------------p2conf.h--------------------------
c
c OTHER GLOBAL CONSTRAINTS ON DIMENSION
c
integer MAXCON, MAXCNT, MAXNB, MAXPAR
integer MAXINF, MAXREG, MAXREC, MAXPTORI, MXBDEP
c.. maximum number of contact nodes
parameter (MAXCON = 500)
c.. maximum number of contacts
parameter (MAXCNT = 10)
c.. maximum number of triangles a point can belong to
parameter (MAXNB = 20)
c.. maximum number material parameters
parameter (MAXPAR = 50)
c.. maximum number of interfaces
parameter (MAXINF = 10)
c.. maximum number of regions
parameter (MAXREG = 1000)
c.. maximum elements in 1-D for temporary tensor-product
c grid info
parameter (MAXREC = 1000)
c.. maximum depth of bias partition
parameter (MXBDEP = 10)
c.. maximum original tensor-product mesh points before
c elimination
parameter (MAXPTORI = 10*MAXPT)
---------------------------------------------------------
In this chapter, four simulation examples will be presented to demonstrate the new capabilities
available in PISCES-2ET. The first two are for the application of carrier temperature analysis. These
are simulations for the substrate current in MOSFET and output characteristics of SOI structures. Both
simulated and measured data are compared. The last two are for the simulation of heterostructures.
Among them the first one shows the electrical simulation of a GaAs/AlGaAs-based LED (for light
emitting diode) with cylindrical symmetry. One unique feature of this example is the existence of the
floating layer in the structure, which poses numerical difficulty. The second example for
heterostructure analysis is the simulation of an AlInAs/GaInAs-based MODFET in which the surface
Fermi-level pinning is shown.
It can be clearly seen that the DUET model provides more reasonable results compared to the
experimental data [42], [43] than does DD model. It should be noticed that during the simulation, one
fitting parameter, which effectively changes the energy relaxation time (or equivalently the energy
relaxation length) used in the expression for energy-dependent ionization rate, is used. The similar
calibration approach has been reported in [22] too. The input deck to PISCES-2ET is listed below for
this device structure. Special attention should be paid to the use of symbolic and model cards. In
the symbolic card parameter engy.ele indicates that during the simulation the electron energy
balance equation should also be solved together with the Poisson’s and carrier continuity equations. In
the model card, three parameters are relevant to the energy-dependent impact ionization (II) model:
imp.tn (but not with imp.tp) specifies that only electron energy dependent II model is used and
for holes the impact ionization rate is still determined by the local field (a default). imp.jt indicates
for energy-dependent II model, the current density rather than the saturation velocity is used to
determine the II rate. Finally, impjt.rat is used to specify the ratio of Lw to τ w v sat (refer to
Eq. (3.39)). The simulation results are shown in Figure 6.2.
title N-MOS for Substrate Current Simulation
$
$ A 0.8 micron N-MOSFET example, Electron-ET model
$
options plotdev=xterm x.scr=3.8 y.scr=3.8
$ mesh
$============================
mesh rect nx=40 ny=28
x.m n=1 l=0 r=1
x.m n=5 l=0.5 r=0.6
x.m n=35 l=1.5 r=1.0
x.m n=40 l=2.0 r=1.7
y.m n=1 l=-0.0152 r=1
y.m n=3 l=0.0 r=1
y.m n=5 l=0.002 r=1
y.m n=9 l=0.010 r=1
y.m n=16 l=0.1 r=1
y.m n=20 l=0.2 r=1
y.m n=22 l=0.3 r=1
y.m n=23 l=0.4 r=1
y.m n=25 l=1.1 r=1
y.m n=28 l=2.0 r=1
$ region and electrodes
$============================
region num=1 ix.l=1 ix.h=40 iy.l=1 iy.h=3 oxide
region num=2 ix.l=1 ix.h=40 iy.l=3 iy.h=28 silicon
$ contact
$============================
contact all neutral
contact num=2 workfunction=4.17
$ initial solution
$============================
symbolic newton carrier=0
method itlim=50 p.tol=1e-5 c.tol=1e-5
model srh fldmob conmob
solve ini
$ ramping vg and vd
$============================
symbolic newton carrier=2 engy.ele
model srh fldmob conmob imp.tn imp.jt impjt.rat=0.69
method itlim=20 p.tol=1e-5 c.tol=1e-5 trap
log ivfile=08_Tn.iv
solve v2=0.2 vstep=0.2 nstep=3 electr=2 proj
solve v2=1.0 outfile=08_Tn.g1 proj
solve v3=0.1 vstep=0.1 nstep=4 electr=3 proj
solve v3=0.6 vstep=0.2 nstep=12 electr=3 proj
solve v3=3.0 outfile=08_Tn.g3 proj
end
Figure 6.1 Input file for simulation of the substrate current for MOSFETs.
-1
10
-2
10
-3
10
I sub/I d
-4
10
MIT
-5 Berkerley
10
0.8 µm
2 µm
Figure 6.2 Simulation results for the substrate current in MOSFET with
two different channel length (2 and 0.8 µ m) and the
comparison is made for 0.8 µ m case between the ET-simulated
and measured results [42], [43]. The upper curves are
simulated using DD model while the lower curves are obtained
plot.1d dop log abs a.x=0 a.y=0.52 b.x=1.48 b.y=0.52 min=14 max=21
$
models conmob intelmob print
mobility ec.crit=8e4 char.int=0.04
material region=2 vsaturation=1e7
symb newton carriers=0
$
method xnorm itlimit=30
$
solve ini
symb newton carriers=1 engy.ele
method trap itlimit=30
$
solve v1=-3 v4=-15 v3=2 vstep=0.5 nstep=10 electr=4
solve v1=-3 v4=64 v3=0.05
end
The partial input deck for SOI with 0.4µm channel length is listed below. Only structural part
is included.
title SOI BERKELEY
$
$ Effective Channel length =.47um
$ 0.47um of silicon on 0.4um oxide substrate
$
mesh width=9.5 rect nx=43 ny=35
$
x.m n=1 l=0.00 r=1.00
x.m n=7 l=0.5 r=0.9
X.m n=12 l=0.6 r=.9
x.m n=22 l=0.915 r=1.15
x.m n=32 l=1.23 r=0.85
x.m n=37 l=1.33 r=1.1
x.m n=43 l=1.83 r=1.1
$
y.m n=1 l=0.0 r=1.0
y.m n=6 l=0.4 r=0.8
y.m n=21 l=0.46 r=1.15
y.m n=33 l=0.52 r=0.8
y.m n=35 l=0.54 r=1.0
$
region num=1 ix.l=1 ix.h=43 iy.l=1 iy.h=6 oxide
region num=2 ix.l=1 ix.h=43 iy.l=6 iy.h=33 silicon
region num=3 ix.l=1 ix.h=43 iy.l=33 iy.h=35 oxide
$
$*********** define the electrodes ************
$ #1-GATE #2-SOURCE #3-DRAIN #4-SUBSTRATE(below oxide)
$
electrod num=1 ix.l=12 ix.h=32 iy.l=35 iy.h=35
electrod num=2 ix.l=1 ix.h=12 iy.l=33 iy.h=33
electrod num=3 ix.l=32 ix.h=43 iy.l=33 iy.h=33
electrod num=4 ix.l=1 ix.h=43 iy.l=1 iy.h=1
$
$*********** define the doping concentrations *****
$
doping uniform conc=6e16 p.type reg=2
doping gauss conc=2e20 n.type characteristic=0.1 start=0.46
+ x.right=0.6 ratio.lateral=0.25 region=2
doping gauss conc=2e20 n.type characteristic=0.1 start=0.46
+ x.left=1.23 ratio.lateral=0.25 region=2
$
interfac qf=20e10
contact num=1 n.poly
...
end
The simulation results are compared with the measured data provided by Hu’s group at UC
Berkeley [44] in Figure 6.4. The agreement between ET model and experiment is good for devices with
both channel lengths. While DD model gives good fit to the measured data for device with 0.47 µ
channel length, it grossly underestimates the drain current for device with 0.12 µ , an indication that
DD models fails to predict the velocity overshoot phenomenon.
guarantees the n-layer is indeed “floating” as far as the potential of the layer is concerned. The effect
of the extra contact on the device characteristics is minor because the blocking layer is doped heavily.
Following are two input files for the simulation. In order to put the contact to the floating layer from
the exterior, the structure is shifted left in the lateral direction such that the rotation axis falls at x = 0 .
For the initial solution, the voltage boundary condition is used first and then the zero-current boundary
condition is switched on. It is observed that before the diode is turned on, i.e., for the bias below 1V,
the numerical (discretization) noise may exceed the magnitude of the simulated current and so the
current conservation law (the sum of the signed terminal currents should be zero) is not well observed.
Because the device behavior before turning on is not a primary concern, bigger bias steps are used and
especially the diode bias is increased directly from 0.7 to 1V. Also note that the n-layer is on the top
of the structure, so the negative cathode bias is incremented. Figure 6.5 lists the first input file for the
simulation for obtaining the initial solution with the voltage boundary condition.
0.006
0.004
Leff=0.12µ
IDS(A)
0.002
Leff=0.47µ measured
DD Model
ET Model
0.000
0.0 0.5 1.0 1.5 2.0
VDS(V)
Figure 6.4 Simulated and measured data [44] for SOI structure
with different channel length.
Figure 6.5 First input deck for the simulation of LED, used to find the initial solution for
voltage boundary condition.
The plots resulted from the running of the above file are shown in Figure 6.6 and Figure 6.7.
Figure 6.6 shows the structure of the device and Figure 6.7 shows the mesh used in the simulation.
Once the initial solution is found with the voltage boundary condition, the current boundary condition
is switched on for the contact to the floating layer in order to ensure it being “floating” (no current
drawn at this terminal). This is implemented in the second file as listed in Figure 6.8.
Figure 6.8 Second input file for LED simulation. The initial solution is obtained from the
voltage boundary condition and is switched to the current boundary condition for
the floating layer contact. The bias step jumps from – 0.7 to – 1 V to avoid the
numerical instability.
Figure 6.9 Simulated I – V characteristics of LED. Notice that the bias is applied to the
top of the structure (n-layer), so it is negative and the current computed below
0.7V is due to the numerical noise.
The numerical convergence behavior during simulation is excellent. The simulated I-V characteristics
is shown in Figure 6.9. The bias ramp is applied to the cathode of the diode, so the magnitude of the
bias increase is in the opposite direction to the x-axis.
$
$ Electrodes: 1 source, 2 gate, and 3 drain
$
elec num=1 ix.l=1 ix.h=4 iy.l=1 iy.h=1
elec num=2 ix.l=20 ix.h=32 iy.l=4 iy.h=4
elec num=3 ix.l=52 ix.h=55 iy.l=1 iy.h=1
$
$ Doping specification
$
doping region=2 unif conc=2e19 n.type
doping region=3 unif conc=2e19 n.type
doping region=4 unif conc=1e15 n.type
doping region=5 unif conc=6e18 n.type
doping region=6 unif conc=1e15 n.type
doping region=7 unif conc=1e15 n.type
doping region=8 unif conc=1e15 n.type
$
$ Background doping
$
doping unif x.l=0.0 x.r=4.0 y.t=0.0 y.b=0.2 conc=5e14 p.type
$
$ Interface trapping for the pinning of Fermi-level
$
deepimp unif x.l=1.625 x.r=1.925 y.t=0.03 y.b=0.031 accep conc=1.0e19
+ eion=0.7
deepimp unif x.l=2.075 x.r=2.375 y.t=0.03 y.b=0.031 accep conc=1.0e19
+ eion=0.7
$
$ Show regions
$
plot.2d bound pause
$
$ Doping and x-mole profiles under the gate
$
plot.1d dop x.s=2.0 x.e=2.0 y.s=0.0 y.e=0.2 log pause
plot.1d xmol x.s=2.0 x.e=2.0 y.s=0.0 y.e=0.2 min=0.45 max=0.49 pause
$
contact num=2 workf=4.8 surf.rec
models temp=300 conmob consrh fldmob auger incompl
$
$ Solving Poisson’s equation only for the equilibrium solution
$
symb newton carriers=0
method itlimit=50 biaspart
$
solve ini
$
$ Switch to full set of semiconductor equations
$
symb newton carr=2
$
Figure 6.10 Input file for the simulation of AlInAs/GaInAs MODFET. The surface Fermi-
level pinning is modeled by specifying the high density of deep level impurities
at the free surface.
The simulation results are shown in the following figures. Figure 6.12 shows the band
diagram at the thermal equilibrium for the Fermi-level pinning due to the high density of the surface
traps and the situation without surface pinning. Finally the simulated drain current vs. the drain voltage
at V GS = – 0.8 V is shown in Figure 6.13.
Figure 6.12 Band diagram along the line segment located at the spacer between the source
and gate in the direction perpendicular to the surface of the MODFET. The
effect of Fermi-level pinning due to the surface traps is shown by comparison
to the free surface.
VGS=-0.8V
2
f 0 ( r, k ) = ---------------------------------------------------------- (A.1)
E ( k ) – E Fn ( r )
1 + exp ----------------------------------
kBT n( r)
We need to transform the function to have the dependence on n, Tn, and ε (the kinetic energy). This
can be done in the following procedure. First with Eq. (A.1), one can show that from the definition of
n (the carrier density) that
E Fn ( r ) – E C ( r )
n = N C ( T n, T L )F 1 ⁄ 2 ------------------------------------- (A.2)
kBT n( r)
where NC has the same expression as in section 2.1 and F 1 ⁄ 2 is the Fermi integral of order one half.
The definition of Fermi integral as used in this manual will be given in APPENDIX B together with
some other related properties. To have more concise notation, we define two symbols:
E Fn – E C
η n = --------------------- (A.3)
kBT n
which is a function of n and Tn as can be seen from Eq. (A.2), and
F 1 ⁄ 2 ( ηn )
γ n ( η n ) = ---------------------
- (A.4)
exp ( η n )
Note that γ n is actually a measure of the carrier degeneracy, i.e., it is always smaller than unity but
approaches the unity for non-degenerate case ( η n → – ∞ ). There are similar expressions for holes, but
η p is expressed as
E V – E Fp
η p = --------------------- (A.5)
kBT p
We can now express the separation of the quasi-Fermi level and band edge using the carrier
concentration as follows:
E Fn – E C n
exp --------------------- = ----------------------------------------- (A.6)
B n k T N C ( T n n ( n, T n )
)γ
Note that throughout this document, we use n and Tn as the fundamental variables. Substituting this
expression in Eq. (A.1), we obtain an equivalent form of the Fermi-Dirac distribution as follows:
2
f 0 ( n, T n, ε ) = ------------------------------------------------------------------------------ (A.7)
N C ( T n )γ n ( n, T n ) ε
1 + -----------------------------------------exp ------------
n k B T n
h h ∂f0
f ( r, k ) = f 0 ( n ( r ), T n ( r ), ε ( k ) ) – τ ( r, k ) ------------*- k ⋅ ∇ f 0 – q ------------*- -------- k ⋅ E n (A.8)
2πm n 2πm n ∂ε
We will further make several assumptions to simplify the above expression. First we assume the
effective mass is constant, i.e. for parabolic band structure. Then the carrier velocity, v , is related to
the wave vector, k , by the following expression:,
h
v = ------------*-k (A.9)
2πm
and the kinetic energy can be written as
* 2 2
h
ε = m v
------------ = -----------------------
2 *
-k
2
(A.10)
2 2 ( 2π ) m
Furthermore, in heterostructures the electric field might be different for electrons and holes depending
on the gradient of the respective band edge. For electrons, we have
1
E n = --- ∇E C (A.11)
q
We assume the relaxation time ( k -dependent τ ) can be simplified to depend on the kinetic energy
only, i.e., τ ( r, k ) = τ ( r, ε ) . Eq. (A.8) can then be written as
∂f0
f ( r, v ) = f 0 ( n, T n, ε ) – τ ( r, ε ) v ⋅ ∇ f 0 – q --------v ⋅ E n (A.12)
∂ε
Now we are ready to evaluate the carrier current and its relevant coefficients. By doing so,
one may also find the exact definition for some commonly used physical parameters such as the
mobility and diffusivity. The carrier current can be evaluated from the distribution function as follows:
* 3
1 mn
∫ ∫ ∫
3
j n = – q vf ( r, k ) -------------dk
3 x dk y dk z = – q v f ( r, v ) ------ dv x dv y dv z = – q v f ( r, v )d v (A.13)
( 2π ) h
3 * 3 –3
where we have used symbol d v = ( m ⁄ h ) dv x dv y dv z and the introduction of ( 2π ) in the integral
in k -space is due to the quantum mechanics consideration. It would be desirable to perform the above
integral over ε if the integrand is the function of r and ε only. This can be done through the following
transformation:
* 3 * 3⁄2
m 2π ( 2m n )
-----n- dv x dv y dv z = ------------------------------
1⁄2 1 –3 ⁄ 2 1⁄2
h ∫ h
3
- ∫ε d ε = ------- ( k B T n )
π
NC ∫ε dε (A.14)
It should be noticed that the first part of the right-hand-side (RHS) in Eq. (A.12) is even in k -
space and the second part is odd in k -space, thus in the integral of Eq. (A.13) the first part of the
distribution function will vanish (by multiplying v the even part becomes odd), and only the second
(odd) part is used in the integration. The result is as follows:
∂f0
j n = q τ ( r, ε )vv ⋅ ∇ f 0 d v – q τ ( r, ε ) --------vvd v ⋅ E n
∫ ∫
3 2 3
(A.15)
∂ε
Note that vv should be understood as the tensor. It is easy to identify the second term on the RHS
represents the drift component because it is proportional to the electric field. The first term is, however,
not clear at this stage and we need to further find out the form of ∇ f 0. We now proceed to find its
expression.
∂f0 ∂f0 ∂ f 0 λn 3 λn ε
∇ f 0 = --------∇n + --------- ∇T n = k B T n -------- – -----∇n + --- ------ – -----------2- ∇T n (A.16)
∂n ∂T n ∂ε n n k T
2 T
B n
Thus Eq. (A.15) becomes
1 3 ∂f0 1 3 ∂f0
j n = q – k B T n λ n --- d vτ ( r, ε ) --------vv ⋅ ∇n + qn – q --- d vτ ( r, ε ) --------vv ⋅ E n
∫ ∫
n ∂ε n ∂ε
(A.17)
31 1 3 ∂f0 11 3 ∂f0
+ qn – --- ------ – k B T n λ n --- d vτ ( r, ε ) --------vv – ------ --- d vτ ( r, ε ) ε --------vv ⋅ ∇T n
∫ ∫
2Tn n ∂ε Tnn ∂ε
This seemly formidably complex expression can greatly be simplified by introducing some familiar
coefficients through identification of origin of each term.
8π 3⁄2
∫ d vvv ∫ε d ε Î
3 *
= -------3- 2m n (A.18)
3h
where Î is the identity matrix. By realizing the relationship of Î ⋅ ∇ = ∇ and by comparing to the
standard drift-diffusion expression of the carrier current, one can write the following compact form:
T
j n = qD n ∇n + qnµ n E n + qnD n ∇T n (A.19)
T
where coefficients D , µ , and D are named as diffusion constant (or diffusivity), mobility, and
thermal diffusivity, respectively. Comparing Eq. (A.19) with Eq. (A.17), one obtains
1 8π ∂f0 3⁄2
µ n = – q --- -------3- 2m n τ ( r, ε ) -------- ε d ε
∫
*
(A.20)
n 3h ∂ε
kBT n
D n = ------------ λ n µ n (A.21)
q
1 3 8π *1 ∂f0 5⁄2
D n = – ------ --- D n + -------3- 2m n --- τ ( r, ε ) -------- ε d ε
∫
T
(A.22)
Tn 2 3h n ∂ε
where Î is the identity matrix. It becomes immediately clear from Eq. (A.20) and Eq. (A.21) that we
have arrived a generalized Einstein relationship between the mobility, µ , and diffusivity, D as
follows:
ε
f 0 ( n, T n, ε ) = -------exp – ------------
n
(A.24)
NC k B T n
there exists a special relationship that
∂f ∂ ∂f0
--------0- = n --------- -------- (A.25)
∂T n ∂T n ∂n
Then we obtain
T ∂D
D n = --------n- (A.26)
∂T n
Note that this simple relationship does not hold for general Fermi-Dirac statistics.
Now we consider the specific form of the relaxation time, τ ( r, ε ) . For the current PISCES-
2ET implementation, we assume the following power dependence:
–ν
τ ( r, ε ) = α ( r, T L ) ε (A.27)
where ν ≥ 0 and α is the remaining part of τ which does not depend on the carrier kinetic energy. To
be more specific, we have explicitly expressed the lattice temperature (T L ) dependence. Also note that
α doesn’t have the units of time.
With the above power dependence of the relaxation time on energy, we can find more direct
T
link between D and µ . Notice that
∞ ∞
∂f
d f 0 = ( 3 – 2ν )Γ --- – ν F 1 ⁄ 2 – ν ( η n ) ( k B T n )
–ν 0 3⁄2 3⁄2–ν 3 3⁄2–ν
∫ ε -------
∂ε
- ε dε = ∫ ε
2
(A.28)
0 0
where integral by parts has been used during the derivation and ν must be no bigger than 3/2 to insure
that no singularity occurs. For the same reasoning, one obtains
∞
–ν ∂ f 0 5 ⁄ 2
d ε = ( 5 – 2ν )Γ --- – ν F 3 ⁄ 2 – ν ( η n ) ( k B T n )
5 5⁄2–ν
∫ε -------- ε
∂ε 2
(A.29)
0
T
D 1 3 kBT n kBT n F3 ⁄ 2 – ν
------n- = – ------ --- ------------ λ n – --- – ν ------------ -----------------
5
(A.30)
µn Tn 2 q 2 q F1 ⁄ 2 – ν
or
T k
D n = --- 5--- – ν F 3⁄2–ν 3
- – --- λ µ
---------------- (A.31)
q 2 F1 ⁄ 2 – ν 2 n n
We now consider a more specific case, that is for the case where the acoustic phonon
scattering is dominant, then p = 1 according to Stratton [1]. Furthermore, for nondegenerate cases,
using property Eq. (B.5), the ratio of Fermi integrals with different orders all become the unity, hence
T
D n even becomes zero. That is to say, the gradient of the carrier temperature doesn’t contribute to the
current flow.
A.3 Heterostructures
For heterostructures, the effective mass would be a function of the position also. In evaluating
∇ f 0, therefore, one needs to take into consideration the variation of the effective mass. This can be
done by adding an extra term in Eq. (A.16) as follows:
1
E p = --- ∇E V (A.33)
q
But for the entire semiconductor region, there is one variable for electrostatic potential, ψ .
So it is desirable to link both E C and E V to ψ . We define for heterostructure the electrostatic potential
as
E vac
ψ = – ---------
- (A.34)
q
where E vac is the energy level for vacuum. The advantage of defining the potential this way is to
warrant it is continuous and differentiable for the electric field has to be finite in magnitude. This is
extremely important as in the heterostructure, the conventional definition of the potential as the
intrinsic Fermi level
EC + EV N C(T L)
- – kT L ln -------------------
E Fi = ------------------- (A.35)
2 N V (T L)
where T L is the lattice temperature and the densities of states for conduction and valence bands are
evaluated at the lattice temperature, is not necessarily continuous. We can then relate the band edges
to the potential using band structure parameters: electron affinity χ and bandgap E g as follows:
E C = – qψ – χ (A.36)
E V = – qψ – χ – E g (A.37)
We are now ready to rewrite the current expressions for electrons and holes. Based on
Eq. (A.32) and Eq. (A.15), one obtains for electrons:
3 kBT n
j n = qD n ∇n + qnµ n F n – --- λ n ------------ ∇lnm n + qnD n ∇T n
* T
2 q
(A.38)
– qnµ n ∇ψ – nµ n ∇χ + --- λ n k B T n ∇lnm n
T 3 *
= qD n ∇n + qnD n ∇T n
2
For holes,
T 3 *
j p = – qD p ∇p – qpD p ∇T p – qpµ p ∇ψ – pµ p ∇ ( χ + E g ) – --- λ p k T p ∇lnm p (A.39)
2 B
The above expressions can greatly be simplified if the quasi-Fermi level instead of the carrier
concentration is used as the driving force for the carrier transport. By introducing the kinetic energy
ε ( k ) = E ( k ) – E C ( r ) in Eq. (A.2), the Fermi-Dirac distribution function can be written as
2 2
f 0 ( r, k ) = ---------------------------------------------------------- = ------------------------------------------------------------------------------ (A.40)
E ( k ) – E Fn ( r ) ε ( k ) + E C ( r ) – E Fn ( r )
1 + exp ---------------------------------- 1 + exp -------------------------------------------------------
k B nT ( r ) kBT n( r)
Thus
∂ f 0 ε ( k ) + E C ( r ) – E Fn ( r )
∇ f 0 ( r, k ) = -------- ∇E C – ∇ E Fn – ------------------------------------------------------- ∇T n (A.41)
∂ε T n
Substituting the above expression into Eq. (A.15) and using Eq. (A.11), one obtains
j n = nµ n ∇E Fn + qnQ n ∇T n (A.42)
where the thermopower Q is defined as
kB 3
Q n = D n + ----- --- λ n – η n µ n
T
(A.43)
q 2
Note that the coefficient for the current component explicitly due to the carrier temperature gradient
in Eq. (A.42) is different from that in Eq. (A.19). But for the isothermal case ( T n = const ), both
expressions are reduced to the standard drift-diffusion (DD) form.
ν
1 x
Γ( ν + 1) 1 + ex – η ∫
F ν ( η ) = --------------------- --------------------
- dx (B.1)
∞
ν – 1 –x
Γ(ν) = ∫x e dx (B.2)
0
Γ --- =
1
π (B.3)
2
and
η
lim F ν ( η ) = e regardless value of ν (B.5)
η → –∞
and
d
------F ν ( η ) = F ν – 1 ( η ) (B.6)
dη
From the above properties, one can derive two useful derivatives:
∂ λ
-----η n ( n, T n ) = ----n- (B.7)
∂n n
where λ n is the symbol of λ ( η n ) with λ defined as
F1 ⁄ 2(η)
λ = ---------------------
- (B.8)
F –1 ⁄ 2 ( η )
∂ γ
-----γ n ( n, T n ) = ----n- ( 1 – λ n ) (B.9)
∂n n
For holes there are completely analogous formulas and are listed below:
∂ λ
------ η p ( p, T p ) = -----p (B.10)
∂p p
∂ γ
------ γ p ( p, T p ) = ----p- ( 1 – λ p ) (B.11)
∂p p
Now we give the expressions for the derivative of η and γ w.r.t. the carrier temperature T .
∂ 31
--------- η n ( n, T n ) = – --- ------ λ n (B.12)
∂T n 2Tn
∂ 31
--------- η p ( p, T p ) = – --- ------ λ p (B.13)
∂T p 2T p
∂ 3 γ
--------- γ n ( n, T n ) = – --- ( 1 – λ n ) -----n- (B.14)
∂T n 2 Tn
∂ 3 γ
--------- γ p ( p, T p ) = – --- ( 1 – λ p ) -----p- (B.15)
∂T p 2 Tp
Considering the limited access to the manuals and reports for the previous versions of PISCES II (i.e.
PISCES-II User’s Manual (1984) [3] and PISCES-IIB Supplementary report (1985) [4]), we list in this
Appendix those equations and expressions referred to in the User’s Manual but not discussed in this
report. Users may find the similar ones in the previous PISCES II reports.
Surface recombination velocities for electrons and holes, v sn and v sp , are given by
** 2
An T
v sn = --------------
qN C
(C.1)
** 2
Ap T
v sp = --------------
qN V
** **
where A n and A p are the effective Richardson constants (p. 389 in [18]) for electrons and holes.
+ ND
N D = --------------------------------------------
( E Fn – E D ) ⁄ kT
-
1 + gD e
(C.2)
- NA
NA = --------------------------------------------
( E A – E Fp ) ⁄ kT
-
1 + gAe
where E D and E A are impurity energies, g A and g D are degeneracy factors for donors and acceptors,
respectively.
Surface mobility reduction factor, gsurf. To model the low-field mobility along the oxide-
semiconductor interface in a simple way, the following formula is used:
q 1⁄2
∆φ b = ----------- E + αE (C.4)
4πε s
where E is the magnitude of the electric field at the interface. Both lowering mechanisms, image force
(the first term on RHS of Eq. (C.4)) and static dipole layer (second term), are considered in the above
model.
[3] M. R. Pinto, C. S. Rafferty and R. W. Dutton, “PISCES-II - Poisson and Continuity Equation
Solver,'' Stanford Electronics Laboratory Technical Report, Stanford University, September
1984.
[5] N. D. Arora, J. R. Hauser, and D. J. Roulston, “Electron and hole mobilities in silicon as a
function of concentration and temperature,” IEEE Trans. Elec. Dev., Vol. ED-29, pp. 292-295,
Feb. 1982.
[6] Z. Yu et al., “Development of doping and temperature dependent mobility model in GaAs using
a robust optimizer,” Private communication, July 1993.
[7] J. M. Dorkel and Ph. Leturcq, “Carrier mobility in silicon semi-empirically related to
temperature, doping and injection level,” Solid-State Elec., Vol. 24, no. 9, pp. 821-825, 1981.
[8] C. Lombardi, S. Manzini, A. Saporito, and M. Vanzi, “A physically based mobility model for
numerical simulation of non-planar devices,” IEEE CAD, Vol. 7, no. 11, pp. 1164-1171, Nov.
1988.
[9] J. T. Watt, Modeling The Performance of Liquid-Nitrogen Cooled CMOS VLSI, Ph.D.
dissertation, Stanford University, May 1989.
[10] A. G. Sabnis and J. T. Clemens, “Characterization of the electron mobility in the inverted <100>
Si surface,” IEDM Technical Digest, pp. 18-21, 1979.
[11] H. Shin, A. F. Tasch, Jr., C. M. Maziar, and S. K. Banerjee, “A new approach to verify and
derive a transverse field-dependent mobility model for electrons in MOS inversion layers,”
IEEE Trans. Elec. Dev., Vol. 36, no. 6, pp. 1117-1124, June 1989.
[12] S. Schwarz and S. Russek, “Semi-empirical equations for electron velocity in silicon: Part II--
MOS inversion layer,” IEEE Trans. Elec. Dev., vol. ED-30, pp. 1634-1639, Dec. 1983.
[14] D. M. Caughey and R. E. Thomas, “Carrier mobilities in silicon empirically related to doping
and field,” Proc. IEEE, pp. 2192-2193, Dec. 1967.
[16] H. R. Yeager, “Circuit-simulation models for the high electron-mobility transistor,” Ph. D.
Thesis, Stanford University, April 1989.
[18] S. M. Sze, Physics of Semiconductor Devices, 1st ed., p. 62, John Wiley & Sons, New York,
1969.
[19] K. Hui, C. Hu, P. George, and P. K. Ko, “Impact ionization in GaAs MESFET’s”, IEEE Elec.
Dev. Lett., vol. 11, no. 3, p. 113, Feb. 1991.
[20] D. Ritter, R. A. Hamm, A. Feygenson, and M. B. Panish, “Anomalous electric field and
temperature dependence of collector multiplication in InP/Ga0.47In0.53As heterojunction
bipolar transistors,” Appl. Phys. Lett. 60 (25), p. 3150, 22 June 1992.
[21] I. Watanabe, T. Torikai, K. Makita, K. Fukushima, and T. Uji, “Impact ionization rates in (100)
Al0.48In0.52As,” IEEE Elec. Dev. Lett., vol. 11, no. 10, p. 437, Oct. 1990.
[24] S. M. Sze, Physics of Semiconductor Devices, 2nd ed., John Wiley & Sons, New York, 1981.
[25] R. S. Muller and T. I. Kamins, Device Electronics for Integrated Circuits, John Wiley & Sons,
New York, 1977.
[26] Sadao Adachi, “Material parameters of In GaxAsyP1-y and related binaries,” J. Appl. Phys. 53
1-x
(12), Dec. 1982, pp. 8775-8792.
[32] S. W. Corzine, et al., “Optical gain in III-V bulk and quantum well semiconductors,” in
Quantum Well Lasers, ed. P. S. Zory, Jr., Academic Press, 1993.
[33] Z. Yu and R. W. Dutton, “SEDAN III – A generalized electronic material device analysis
program,” Tech. Rep. Stanford University, 1985.
[36] S. E. Laux, “Techniques for small-signal analysis of semiconductor devices,” IEEE Trans. Elec.
Dev., Vol. ED-32, No. 10, pp. 2028-2037, 1985.
[37] Z.-Y. Wang, K.-C. Wu, and R. W. Dutton, “An approach to construct pre-conditioning matrices
for block iteration of linear equations,” IEEE Trans. CAD, Vol. 11, No. 11, pp. 1334-1343,
1992.
[38] R. W. Freund, RIACS Tech. Rep. 91.18, NASA Ames Res. Ctr., Sept. 1991.
[39] R. W. Freund and N. M. Nachtigal, RIACS Tech. Rep. 90.51, NASA Ames Res. Ctr., Dec.
1990.
[40] K. Wu, Z. Yu, L. So, R.W. Dutton, and J. Sato-Iwanaga, “Robust and Efficient AC Analysis of
High-speed Devices,” IEDM Technical Digest, pp. 935-938, San Francisco, Dec. 1992.
[42] T. Y. Chan, P. K. Ko, and C. Hu, “A simple method to characterize substrate current in
MOSFET’s,” IEEE Elec. Dev. Lett., Vol. 5, no. 12, pp. 505-507, Dec. 1984.
[46] D. J. Rose and R. E. Bank, Global Approximate Newton Methods, Numerische Mathematik, 37,
pp. 279-295, 1981.
CARD FORMAT
PISCES-2ET takes its input from a user provided ASCII file. Each line is
a particular statement (or card), identified by the first word on the card
(i.e., card name). The remaining parts of the line are the parameters of that
statement. The words on a line are separated by blanks or tabs. If more
than one line of input is necessary for a particular statement, it may be
continued on subsequent lines by placing a plus sign (+) as the first non-
blank character on the continuation lines. Card and parameter names do
not need to be typed in full; only enough characters to ensure unique
identification is necessary.
In the card descriptions which follow, the letters required to identify a pa-
rameter (i.e. a word) are printed in upper case, and the remainder of the
word in lower case. For each card, parameters are first presented in a SYN-
TAX section and are then described in a PARAMETER section.
CARD SEQUENCE
• The MESH card must precede all other cards, except TITLE,
COMMENT, and OPTIONS cards.
• When defining a rectangular mesh, the order of specification is
as follows:
MESH
X.MESH (all)
Y.MESH (all)
ELIMINATE
SPREAD
REGION
ELECTRODE
ELIMINATE and SPREAD cards are optional but if they occur
they must be in that order.
• DOPING cards must follow directly after the mesh definition.
• Before a solution, a symbolic factorization is necessary. Unless
solving for the equilibrium condition (initial solution), a
previous solution must also be loaded to provide an initial
guess.
• Any CONTACT cards must precede the SYMBOLIC card.
• Physical parameters may not be changed using the MATERIAL,
CONTACT, or MODEL cards after the first SOLVE or LOAD
The execution of the program is line based, i.e. whenever a line is read in it
is processed immediately. This feature allows users to run PISCES interac-
tively.
CONVENTIONS
{ par1 | par2 }
states that either par1 or par2 may be specified but not both.
All file names are currently set to no more than 20 characters in length.
Table M.1 lists the abbreviation of symbols for units in this manual.
Abbr. Units
K Kelvin
J Joule
s second
µ micrometer
A ampere
W watt
C coulomb
COMMAND CHECK
SYNTAX
PARAMETERS
Infile
A character string specifying the name of the solution file to compare.
(Default: null)
Meshfile
A character string as the name of the file containing the mesh for the
solution specified by Infile. (Default: null)
Samemesh
A logical flag indicating that the solution in Infile used the same mesh as
the current solution. (Default: false)
EXAMPLES
COMMAND COMMENT
SYNTAX
PARAMETERS
None.
EXAMPLES
COMMAND CONTACT
SYNTAX
PARAMETERS
ALL
A logical flag to indicate the same properties defined in this card are to be
applied to all electrodes. (Default: false)
ALPha
A real number for the linear, dipole barrier lowering coefficient ( α in
Eq. (C.4)) in units of cm. (Default: 0.0)
ALUminum
A logical flag for using the work function of aluminum (4.17 eV).
(Default: false)
Barrierlow
A logical flag to turn on the barrier lowering mechanism (Eq. (C.4)).
(Default: false)
CApacitance
A real number for the capacitance per unit width (in z-direction) of a
lumped capacitor to be attached to the contact, in units of F µ -1.
(Default: 0.0)
COn.resist
A real number for the distributed contact resistance, ρ c , in units of Ω cm2.
The actual resistance associated with a certain area, A, of the contact is
computed as R = ρ c ⁄ A . (Default: 0.0)
CUrrent
A logical flag to impose the current boundary condition. (Default: false)
MO.disilicide
A logical flag for using the work function of molybdenum silicide
(4.80 eV). (Default: false)
MOLybdenum
A logical flag for using the work function of molybdenum (4.53 eV).
(Default: false)
NEutral
A logical flag to indicate that the work function is to be calculated from
the doping level in the semiconductor at the contact assuming charge
neutrality and thermal equilibrium. Note that neutral stands for charge
neutral. (Default: true)
N.polysilicon
+
A logical flag for using the work function of n silicon material, which is
assumed to be the electron affinity for silicon (4.17 eV). (Default: false)
NUmber
An integer for the electrode sequence number which must previously be
defined in the ELEctrode card. (No default)
Resistance
A real number for the resistance per unit width (in the z-direction) of a
lumped resistor to be attached to the contact, in units of Ω µ .
(Default: 0.0)
Surf.rec
A logical flag to indicate that finite surface recombination velocities are to
be used at the respective contact. (Default: false)
P.polysilicon
+
A logical flag for using the work function of p silicon material, which is
assumed to be the sum of the electron affinity and bandgap of silicon
(4.17 + Egap eV). (Default: false)
TH.resist
A real number for the thermal resistance, as defined in Eq. (2.28), per unit
width (in z-direction) attached to the contact, in units of K s J-1 µ .
(Default: 0.0)
TU.disilicide
A logical flag for using the work function of tungsten silicide (4.80 eV).
(Default: false)
TUNgsten
A logical flag for using the work function of tungsten (4.63 eV).
(Default: false)
VSURFN, VSURFP
Real number parameters for electron and hole surface recombination
velocities, respectively, in units of cm/s. (Defaults: as calculated according
to Eq. (C.1))
Workfunction
A real number for the work function to be used, in units of eV. (Default:
using flag NEutral)
EXAMPLES
COMMAND CONTOUR
SYNTAX
LIne.type = <integer>
ABsolute = <logical> LOgarithm = <logical>
{ X.compon = <logical> | Y.compon = <logical> }
PAuse = <logical> COLor = <logical>
C1.color = <integer> C2.color = <integer>
C3.color = <integer> C4.color = <integer>
C5.color = <integer> C6.color = <integer>
C7.color = <integer> C8.color = <integer>
C9.color = <integer> C0.color = <integer>
PARAMETERS
ABsolute
A logical flag to specify that the absolute value of the plotted quantity be
taken. (Default: false)
ALPHAN, ALPHAP
Logical flags for impact ionization rates due to electrons and holes ( α n
and α p in Eq. (3.37)), respectively. (Defaults: false)
COLor
A logical flag to specify that color fill, as opposed to simple lines, should
be used to delineate contours. (Default: false)
CONduc.b
DEl.value
A real number parameter to specify the interval between values of adjacent
contours. (No default)
DOping
A logical flag for net doping concentration. (Default: false)
E.field
A logical flag for electric field. (Default: false)
ELectrons
A logical flag for electron concentration. (Default: false)
Flowlines
A logical flag for current flow lines. (Default: false)
GEN.Elec, GEN.Hole
Logical flags for impact ionization generation rates due to electrons and
holes (first and second terms on the right hand side of Eq. (3.37)),
respectively. (Defaults: false)
Holes
A logical flag for hole concentration. (Default: false)
Impact
A logical flag for generation rate due to impact ionization (G in Eq.
(3.37)). (Default: false)
LIne.type
An integer parameter to define the type of plot line. (Default: 1)
LOgarithm
A logical flag for using logarithmic instead of original value of the
quantity to plot. For rapidly varying quantities, the logarithmic value is
often more revealing. Since many of the quantities may be negative, the
program actually uses
MIn.value, MAx.value
Real number parameters to specify the minimum and maximum values of
the quantity to be plotted, respectively. If the value is logarithmic, the
minimum and maximum should be given as the logarithmic bounds.
(Defaults: the actual minimum and maximum values of the quantity to be
plotted over the device)
NContours
An integer parameter to specify the number of contours to be plotted. This
is an alternative to the above specification. (No default)
NET.CArrier
A logical flag for net carrier concentration. (Default: false)
NET.CHarge
A logical flag for net charge density. (Default: false)
PAuse
A logical flag to cause the program to stop at the end of the plot so that a
hard copy may be made before continuing. Execution can be resumed by
hitting a carriage return. (Default: false)
POtential
A logical flag for electrostatic potential which is usually defined as the
intrinsic Fermi level. (Default: false)
QFN, QFP
Logical flags for electron and hole quasi-Fermi levels, respectively.
(Defaults: false)
Recomb
A logical flag for net recombination rate. (Default: false)
VAlenc.b
A logical flag for potential corresponding to the valence band edge.
(Default: false)
VELO.Elec, VELO.Hole
Logical flags for electron and hole velocities, which are derived from the
carrier concentration and current density, respectively. (Defaults: false)
X.compon, Y.compon
Logical flags for taking x and y components of a vector quantity,
respectively. (Default: false)
EXAMPLES
COMMAND DEEPIMPURITY
This card has essentially the same features as the DOping card except it
allows users to specify the impurity ionization energy as well. Its main
application is for trap DC analysis as shown in the example for MODFET
simulation.
SYNTAX
Eioniz = <real>
PARAMETERS
Eioniz
Real number parameter to specify the ionization energy for (deep)
impurity level, in units of eV. It is EC - ED for donors and EA - EV for
holes. (Default: same as EDb and EAb in MAterial card)
EXAMPLES
COMMAND DOPING
The DOPING card specifies the impurity doping in selected regions of the
device.
SYNTAX
REgion = <integer>
profile specification:
<profile>
COncentrat = <real> JUnction = <real>
J.Conc = <real> SLice.lat = <real>
or
COncentrat = <real>
{ N.type or DONor = <logical> |
P.type or ACceptor = <logical> }
If <profile type> SUprem3 or OLD.Suprem3:
<Input file>
Infile =<filename>
<dopant>
Boron = <logical> PHosphor = <logical>
ARsenic = <logical> ANtimony = <logical>
The selected dopant file will be extracted from the
SUPREM-III save file
<Two-dimensional spread>
DIrection = <x or y> STart = <real>
RAtio.lat = <real>
<Input file>
Infile =<filename>
<dopant type>
N.type or DONor = <logical>
P.type or ACceptor = <logical>
The ASCII concentrations in Infile are read and added
to (N.type) or subtracted from (P.type) the impurity
profiles.
<Two-dimensional spread>
DIrection = <x or y> STart = <real>
RAtio.lat = <real>
If <profile type> = S4geom or SImpldop:
Infile = <filename>
save:
Outfile = <filename>
PARAMETERS
ACceptor (P.type)
A logical flag to specify p-type doping impurity (i.e. acceptors).
(Default: false)
AScii
A logical flag for input doping file being in ASCII. If AScii is specified
with SUprem3, then an ASCII SUPREM-III export file is expected.
AScii without SUprem3 allows for input a simple ASCII data file
containing information of concentration versus depth. The format of the
ASCII input file is a depth in µ followed by a concentration in cm-3-- one
pair per line. By convention, positive concentrations refer to donors (n-
type), while negative concentration values refer to acceptors (p-type).
(Default: false)
CHaracter
A real number parameter for the principal characteristic length used in
computing profiles specified by either Erfc or Gaussian, in units of µ.
(No default)
COncentrat
A real number parameter for the peak doping concentration, in units of
cm-3. For a Uniform profile, it is the value of the doping level. (No
default)
DIrection
A character of either x or y. Along with STart and RAtio.lat, it specifies
where to locate a one-dimensional profile in the two-dimensional device
and how to extend it to the second dimension. DIrection is the axis along
which the profile is to be directed. (Default: y)
DONor (N.type)
A logical flag to specify n-type doping impurity (i.e. donors).
(Default: false)
DOSe
A real number parameter for the total dose, in units of cm-2. (No default)
ERFC.lat
A logical flag to specify ERFC doping distribution in the lateral direction.
(Default: false)
Infile
A character string for the name of the input file from which the doping is
either directly obtained or interpolated onto an existing mesh. If AScii is
also specified, the doping profile specified in Infile is added to the
previous, if any, impurity profile. (Default: null)Sc
J.conc
A real number parameter for the concentration at the junction, in units of
cm-3. (Default: COncentrat / 100)
JUnction
A real number parameter for the location of the junction, in units of µ. it
must be located in semiconductor region, outside the constant doping box.
When JUnction is used, the program computes the characteristic length
by examining the doping at a point half way between the ends of the
constant box and at the given depth. If some other lateral position is
desired for the computation, use the parameter SLice.lat. (No default)
Lat.char
A real number for the characteristic length of the distribution in the lateral
direction as opposed to that in the principal direction (CHaracter). (No
default)
OLd.Suprem3
A logical flag to read the doping profile from the earlier version of
SUPREM-III process simulation program. A binary structure file
specified by Infile is read in. The defaults for the constant doping box are
set up as a line, parallel to the surface, and located at STart.
(Default: false)
OUtfile
A character string for the name of a binary file in which all the DOping
cards in the present file will be saved. The first DOping card should have
the OUtfile parameter, so that the doping information on it and all
subsequent DOping cards are saved in that file. The file can be re-read
after regridding to calculate the doping profile on a new mesh.
(Default: null)
PEak
A real number parameter specifying the peak position of a doping profile,
in units of µ. (Default: 0.)
RAtio.lat
A real number parameter governing the profile outside the constant doping
box. The lateral profile is assumed to have the same form as the principal
one, but is shrunk/expanded by the factor RAtio.lat. (Default: 0.8)
REgion
An integer parameter for the sequence number of the region where doping
profile is to be added. This is an optional parameter. Multiple regions may
SImpldop
A logical flag specifying use of the doping profile from a SIMPL-2 file
(rectangular grid) to be interpolated onto an existing PISCES mesh.
(default: false)
SLice.lat
A real number parameter, in units of µ, to specify the location in the lateral
direction (as opposed to the principal one), at which the junction depth
specified by JUnction is to be searched along the principal direction. (No
default)
STart
A real number parameter to specify where a one-dimensional doping
profile is to be located along the direction specified by DIrection, used
together with RAtio.lat also. (Default: 0.0)
SUprem3
A logical flag to read the doping profile from the “export” file saved during
a process simulation using the late release of SUPREM-III. The default
is to read binary export files. If the AScii parameter is also specified, then
ASCII SUPREM-III export files will be expected. The defaults for the
constant doping box are set up as a line, parallel to the surface, and located
at STart. (Default: false)
S4geom
A logical flag to read the doping profile from the “geometry” file saved
during SUPREM-IV 2D process simulation. (Default: false)
Table M.1 Default bounding parameters for constant doping box. Uniform
refers to the Uniform profile while x- and y-directions apply to one of
Erfc, Gaussian, SUprem3, OLD.Suprem3, or AScii profiles when the
principal direction is specified as such.
Default
Parameter All other profile type
Uniform
x - direction y - direction
X.Left –∞ SP –∞
X.Right ∞ X.Left ∞
Y.Top –∞ –∞ Y.Bot
Y.Bottom ∞ ∞ SP
EXAMPLES
COMMAND ELECTRODE
SYNTAX
Number = <integer>
position:
Thermal = <logical>
PARAMETERS
Contact
A logical flag for placing an electrode on the top of the surface of the
semiconductor region (such as at the insulator-semiconductor interface) in
Number
An integer number for the electrode. There may be up to ten electrodes,
numbered 1, 2, 3, ... , 9, 0. They may be assigned in any order, but if there
are N electrodes, none can have an electrode number above N.
SUBstrate
A logical flag for placing an electrode at the bottom of the semiconductor
region in the device structure. If Y.Low is specified, then the electrode has
thickness of Y.Low extending to the semiconductor region.
(Default: false)
SURface
A logical flag for placing an electrode on the top (i.e., surface) of the
device structure.(Default: false)
SYmmetri
A logical flag to place an electrode which is symmetric with the one
specified in this card. The device must have a symmetric structure
(Default: false)
Thermal
A logical flag to specify that the electrode is a thermal contact and there is
no equal-potential requirement for those nodes on the same “electrode”
EXAMPLES
COMMAND ELIMINATE
SYNTAX
PARAMETERS
X.direction, Y.direction
Logical parameters to determine whether to eliminate grids along vertical
(Y.direction) or horizontal (X.direction) lines. One must be chosen.
EXAMPLES
The END card specifies the end of a set of PISCES input cards. The END
card may be placed anywhere in the input deck; all input lines below the
occurrence of the END card will be ignored. If an END card is not
included, all cards in the input file are processed.
SYNTAX
ENd or Quit
PARAMETERS
None.
EXAMPLES
END
COMMAND EXTRACT
SYNTAX
Outfile = <filename>
PARAMETERS
Contact
An integer as the number of the contact for which electrode quantities, e.g.
current and (metal) charge, are integrated in its portion intersecting the
Metal.charge
A logical flag to indicate that integrated charge on a contact or its portion
is to be calculated. This flag is useful for studies such as capacitance.
(Default: false)
N.Current, P.Current
Logical flags for computing electron and hole currents, respectively,
through an electrode or its portion. (Defaults: false)
N.Resist, P.Resist
Logical flags for the resistance of a cross section due to electrons and
holes, respectively. The resistance of a cross section is defined as the one
for the resistor with unit thickness in the third dimension and using the
cross section as contacts. That is (for electrons),
1
R n = --------------------- (M.2)
q ∫ nµ n ds
Outfile
Character string for the name of an optional ASCII output file to which the
result and bias information are to be written. (Default: null)
Regions
An integer to identify region(s) as defined in Region card. The
integration will be conducted over only those nodes which fall within both
the bounds specified by X.MIn, X.MAx, Y.MIn, Y.MAx and the particular
set of regions specified by Regions. (No default)
EXAMPLES
EXTRACT P.RESIST
Extracts the resistance of a p-type line diffused into a lightly doped n
substrate. Since the p-conductivity of the substrate is negligible, the
bounds of the integration can include the whole device.
COMMAND IMPACT
The IMPACT card specifies the impact ionization model used for (local)
electric-field dependence. For many devices, the impact ionization model
for continuity equations allows the accurate prediction of avalanche
breakdown. The current models are for Si only. See also, relevant
parameters in the MODELS and CONTOUR cards for both field and
(carrier) energy dependent impact ionization models.
SYNTAX
PARAMETERS
Break
A logical flag to calculate the breakdown voltage using the ionization
integral. (Default: false)
Crowell
A logical flag to use the analytical expression proposed by Crowell and
Sze for temperature dependent impact ionization rates [45]. There are four
parameters in this expression: ionization energy, ε i = 1.5E g , Raman
optical phonon energy, ε r , carrier mean free path for optical phonon
generation, λ , and the low temperature limit of the mean free path, λ 0 . λ
for electrons and holes can be altered by users through LAMDAE and
LAMDAH in this card. The other parameters have values as used in [45].
If Crowell is not specified, the model as specified in Eq. (3.38) and
Table 3.11 is used. (Default: false)
Ionpath
A logical flag to print the nodes, electric field, carrier ionization rates
along the integration path for calculation of breakdown voltage.
(Default: false)
LAMDAE, LAMDAH
Real number parameters for electron and hole mean free paths ( λ in [45]),
–7
respectively, in units of cm. (Default: LAMDAE = 6.2 ×10 , LAMDAH =
–7
3.8 ×10 at 300 K [45])
Monte
Character string for the name of the output solution file from a Monte
Carlo (MC) solver. A non-null string serves also as the flag to indicate that
the impact ionization rates ( α ’s) are to be calculated from this MC
solution file. (Default: null)
EXAMPLES
COMMAND INCLUDE
SYNTAX
A character string
PARAMETERS
None.
EXAMPLES
INCLUDE MAT.ini
COMMAND INTERFACE
SYNTAX
Number = <integer>
parameters:
PARAMETERS
Number
An integer for the sequence number of the interface. The maximum
number is currently set to 10. (No default)
Qf
Real number for the interfacial fixed charge density, in units of cm-2.
(Default: 0.0)
S.N, S.P
Real number parameters for electron and hole surface recombination
velocities, respectively, at the semiconductor surface, in units of cm/s.
(Defaults: 0.0)
EXAMPLES
COMMAND LOAD
The LOAD card loads previous solutions from files for plotting or as
initial guess for present solution at other bias point.
SYNTAX
PARAMETERS
Ascii
A logical flag to specify that any files read in or written to by Load should
be ASCII rather than binary. (Default: false)
Difference
A logical flag to indicate that the difference between two solutions
(IN1file and IN2file) is to be analyzed. The difference may be stored by
specifying Outdiff and can only be used for the purpose of data plotting
or extraction. (Default: false)
No.check
A logical flag to prevent PISCES from checking material parameter
difference between the loaded file(s) and those specified in the current
input file. Checking will never be done for loading of ASCII solution
file(s). (Default: false)
Outdiff
Character string for name of file to which the difference from comparison
of loaded solution files (flag Difference) is written to. (Default: null)
Restore
A logical flag to restore the material parameters and models in the loaded
file as the current ones. (Default: false)
EXAMPLES
LOAD INF=SOL.IN
Specifies that a single solution file called SOL.IN should be loaded.
COMMAND LOG
The format of the data saved in AC log file is as follows. The first line in
the file records the number of terminals for the device under simulation.
The lines starting with “*” record the number of the activation terminal,
the magnitude of the AC input signal (sinusoidal), frequency, and DC bias
(terminal voltages if lumped external resistors exist). The lines without
“*” record, in sequence, the number of the activation terminal, AC
conductances, capacitances, and magnitude of admittances.
SYNTAX
PARAMETERS
Acfile
Character string for the name of file to which the AC simulation results are
to be written. (Default: null, meaning no AC data saved)
Ivfile or Outfile
Character string for the name of file to which simulated I-V information is
to written. (Default: null, meaning no DC data saved)
EXAMPLES
COMMAND MATERIAL
SYNTAX
PARAMETERS
AFfinity
Electron affinity, in eV.
ARICHN, ARICHP
** **
Richardson constants for electrons and holes ( A n and A p in Eq. (C.1)),
respectively, in units of A K-1 cm-2. (Default: Table M.2).
AUGN, AUGP
Auger coefficients for electrons and holes ( c n and c p in Eq. (3.47)),
respectively, in units of cm6 s-1. (Default: Table M.2)
EAb, EDb
Acceptor and donor ionization energies (E A – E V and E C – E D in Eq.
(C.2)), respectively, in units of eV. (Default: Table M.2)
EGAlpha
α in Eq. (4.2) of units eV K-1.
EGBeta
β in Eq. (4.2) of units K.
EG300
Energy bandgap at 300 K (Eq. (4.2)), in eV.
ETrap
Trap level = Et -Ei, in units of eV. (Default: Table M.2)
GCb, GVb
Degeneracy factors for the donor and acceptor levels ( g D and g A in
Eq. (C.2)), respectively, unitless. (Default: Table M.2)
G.surface
Surface mobility reduction factor (gsurf in Eq. (C.3)), unitless. (Default:
Table M.2)
KP.300
Thermal conductivity at 300 K ( κ 0 in Eq. (3.36)), in units of W cm-1 K-1.
(Default: Table 3.10)
KP.Alpha
Temperature coefficient for the thermal conductivity ( α in Eq. (3.36)),
unitless. (Default: 1.2)
MUN, MUP
Low-field mobilities for electrons and holes, respectively, in cm2 V -1 s-1.
NC300, NV300
Effective densities of states for conduction and valence bands,
respectively, at 300 K, in cm-3.
NSRHN, NSRHP
Reference concentration in SRH recombination formula (N SRH in
Eq. (3.45)) for electrons and holes, respectively, in cm-3.
(Default: Table M.2)
NUmber or Region
An integer for the sequence number of the region where any of three
material parameters: static dielectric constant (i.e. Permittivity), thermal
conductivity at room temperature (KP.300) and its temperature coefficient
(KP.Alpha) can be changed from region to region. Note that all other
material parameters can only be changed for the entire device (base
material). (No default)
Permittivity
Static dielectric constant, unitless.
TAUN0, TAUP0
Minority carrier lifetimes for electrons and holes as used in SRH formula
( τ 0 in Eq. (3.45)), in s. (Default: Table M.2)
Vsat
Carrier saturation velocity, in cm/s. (Default: see Table 3.9 and Table 4.1)
Insulator Permittivity
Silicon dioxide 3.9
Silicon nitride 7.5
Sapphire 12.0
EXAMPLES
COMMAND MESH
The MESH card either initiates the mesh generation phase or reads a
previously generated mesh.
SYNTAX
<Previous>
Infile = <filename> ASCII.In = <logical>
<Rectangular>
Rectangular = <logical>
NX = <integer> NY = <integer>
Diag.flip = <logical>
<Geometry>
Geometry = <logical> Infile = <filename>
Flip.y = <logical> SCale = <integer>
output files:
smoothing key:
SMooth.key = <integer>
dimension:
PARAMETERS
ASCII.In
A logical flag to indicate that the input mesh file is an ASCII rather than a
binary file. (Default: false)
ASCII.Out
A logical flag to indicate the mesh file, as specified by OUTFile, to be
saved is an ASCII one. Otherwise it will be binary. (Default: false)
Cylindrical
This logical flag specifies that the mesh, whether generated in the current
run or read in from a file, is to be rotated about the right edge of the
rectangular simulation region to permit the simulation of cylindrically
symmetrical devices. This information will not be written to the mesh file;
it must be specified in the input deck for each PISCES run. (Default: false)
Diag.flip
A logical flag to flip the diagonals in a rectangular mesh about the center
line of the mesh. If not set, all the diagonals will be in the same direction.
(Default: false)
Elec.bot
A logical flag to add an electrode to the bottom of the simulation region
when the mesh file read in is a geometry file generated by a 2D process
simulator such as SUPREM-IV. One of the possible applications of this
flag is to add a substrate contact for a MOSFET structure. Currently, the
electrode number assigned to this contact is hard-wired to 4 in the code.
(Default: false)
Flip.y
A logical flag which reverses the sign of the y-coordinate. (Default: false)
Geometry
A logical flag to indicate that the input ASCII mesh file as specified by
Infile was generated by either a 2D process simulator or by an external
mesh generator. (Default: false)
Infile
The name of the previously generated mesh file. If flag ASCII.In is set,
this file is an ASCII file. Otherwise, it is binary. If flag Geometry is set,
this mesh file, which must be an ASCII one, was generated by either a 2D
process simulator such as SUPREM-IV or by an external mesh generator.
(Default: null)
NX, NY
Integer parameters for numbers of grids in the x- and y-directions for a
rectangular mesh, respectively. (No defaults)
OUT.ascii
This logical flag indicates that an ASCII mesh file will be saved for later
editing by an external grid editor. See Appendix B of [3] for details of the
format. (Default: false)
OUTFile
Specifies the name of the mesh file to be saved. The saved mesh file is
either ASCII or binary depending on the flag ASCII.Out and can be read
in by a later PISCES run. (Default: null)
Poly.elec
A logical flag to convert the region with material type of polysilicon to an
electrode. The condition and application scope of this flag are the same as
those of Elec.bot. Because of the potential difference in integer
representation of the polysilicon region for different process simulators,
this flag is not guaranteed to work properly but it has been tested using one
of the commercial versions of SUPREM-IV. Currently the electrode
number assigned to this contact is hard-wired to 3 in the code. (Default:
false)
Rectangular
A logical flag which initiates the generation of a rectangular mesh.
(Default: true)
SCale
An integer factor by which all the coordinates read in are multiplied.
(Default: 1)
SMooth.key
This integer parameter causes mesh smoothing as described in Section 4.6
of [3]. The digits of the integer are read in reverse order and decoded as
follows:
1. Triangle smoothing, maintaining all region boundaries fixed.
2. Triangle smoothing, maintaining only material boundaries.
3. Node averaging.
Options 1 and 3 are the most common; 2 is used only if a device has
several regions of the same material and the border between the different
regions is unimportant. (No default)
Width
This real number parameter specifies the width of the simulation region,
i.e., the dimension in the z-direction, in units of microns. When this
parameter is used, all units for currents become amperes instead of A/µ,
otherwise simulated terminal currents are considered as current densities
per micron in width. (Default: 1.0)
EXAMPLES
COMMAND METHOD
SYNTAX
<Gummel method>
Gemmul iteration)
SInglepois = <logical> ICcg=<logical>
LUcrit or LU1crit = <real> LU2crit = <real>
Maxinner = <integer> ACCElerat = <logical>
ACCSTArt = <real> ACCSTOp = <real>
ACCSTEp = <real>
<Newton method>
PARAMETERS
AUtonr, Nrcriter
These two parameters are for implementing an automated Newton-
Richardson procedure which attempts to reduce the number of LU
decompositions per bias point. The logical flag AUtonr indicates that this
algorithm is to be used. Real number parameter Nrcriter is the ratio by
which the norm from the previous Newton loop must go down in order to
be able to use the same Jacobian (i.e., LU decomposition) for the current
Newton loop. This is strongly recommended for full Newton iteration.
(Defaults: AUtonr, false; Nrcriter, 0.1)
DAMPEd
A logical flag to indicate the use of a more sophisticated damping scheme
proposed by Bank and Rose [46] (this is the recommended option,
particularly for large bias steps). (Default: false)
DAMPLoop
An integer for the maximum number of damping loops allowed to find a
suitable damping coefficient. (Default: 10)
DElta
Real number parameter for the threshold value to determine the damping
factor for ∆ψ, must be between 0 and 1. (Default: 0.5)
DFactor
Real number parameter for the factor which serves to increase the initial
damping coefficient for the next Newton loop. (Default: 10.0)
DT.min
Real number parameter used in transient simulation for the minimum
– 25
time-step allowed, in seconds. (Default: 1 ×10 )
DVlimit
Real number parameter to limit the maximum update in ψ for a single
loop, in units of the thermal voltage kT/q. (Default: 0.1)
Extrapolate
A logical flag used in transient simulation. Uses a second-order
extrapolation to compute initial guess for the successive time-step.
(Default: false)
Fix.qf
A logical flag to fix the quasi-Fermi potential of each non-solved for
carrier to a single value, instead of picking a value based on local bias (see
Section 5 of Chapter 2 [3] and the ‘‘P.bias'' and ‘‘N.bias'' parameters on
the SOlve card). (Default: false)
ICcg
A logical flag to choose whether or not to use iteration to solve the multi-
Poisson loops. It should be set whenever doing multi-Poisson.
(Default: false)
ITlimit
An integer parameter for the maximum number of allowed outer loops
(i.e., Newton loops or Gummel continuity iterations). (Default: 15)
LImit
A logical flag to indicate that the convergence criterion should be ignored,
and iterations are to proceed until ITlimit is reached. (Default: false)
L2norm
A logical flag to indicate use of L2 norm in transient simulation. It
specifies that the error norms be L2 as opposed to infinity norms for
calculating the time-steps. (Default: true)
Maxinner
An integer to set the maximum number of ICCG iterations. (Default: 25)
PRint
A logical flag to print the terminal fluxes/currents after each continuity
iteration; if this parameter is not set, the terminal fluxes/currents are only
printed after the solution converges. (Default: false)
Rhsnorm
A logical flag. If selected, the Poisson error is measured in C µ-1 and the
continuity error in A µ-1. P.toler then defaults to 1 ×10 C µ-1, and
– 26
SInglepois
A logical flag to indicate that only a single Poisson iteration is to be
performed per Gummel loop as opposed to the default where the
TAuto
A logical flag used in transient simulation to force the program to select
time-steps automatically from the local truncation error estimates. Note
that automatic time-stepping is the default for the second-order
discretization but is not allowed for the first-order scheme. (Default: true
when flag 2ndorder set)
TOl.time
Real number parameter used in transient simulation for the maximum
–3
allowed local truncation error, unitless. (Default: 5 ×10 )
Xnorm
A logical flag to indicate that the error norm in Poisson updates is
measured in units of kT/q, and that in carrier updates is measured relative
to the local carrier concentration. In this case the default value for both
–5
P.toler and C.toler is 1 ×10 . (Default: true)
2ndorder
A logical flag used in transient simulation to specify that the second-order
discretization of Bank, et al. [47] be used as opposed to the first-order
backward difference method. (Default: true)
EXAMPLES
and the Poisson error tolerance should be 1 ×10 C µ-1. Note that because
– 30
Xnorm defaults to true, Xnorm must be turned off to use the rhs norm as
a convergence criterion. If Xnorm = false had not been specified, the rhs
norm and the update norm would have both been printed, but only the
update norm would have been used to determine convergence.
COMMAND MOBILITY
The MOBILITY card specifies the values for the parameters in the
mobility models.
SYNTAX
HOle = <logical>
region:
BULk = <logical>
numerical parameters:
PARAMETERS
Alpha
Real number for the coefficient α in Intel’s local field dependent mobility
model of transverse field reduction (Eq. (3.16)). (Defaults: 1.02 for
electrons and 0.95 for holes)
BNL, BPL
Real number parameters for coefficient B in Eq. (3.18) (Lombardi’s
mobility model) for electrons and holes, respectively. (Defaults: BNL,
7 7
4.75 ×10 ; BPL, 9.93 ×10 )
BUlk
A logical flag to indicate that parameters specified in the present card
apply to the bulk as opposed to the surface mobility model. These
parameters include: EXP1.tem, EXP.Conc, EXP2.tem, Ncrit, UMAx,
and UMIn (Default: false, meaning that the surface mobility parameters
are applied)
CHar.int
Real number parameter of λ in Eq. (3.5) to specify the decay length for
the complementary error function used for the transition from bulk to
surface mobility, in units of µ. (Default: 0.03)
CNExp, CPExp
Real number parameters for β in Eq. (3.18) (Lombardi’s mobility model)
for electrons and holes, respectively. (Defaults: CNExp, 0.125 ;
CPExp, 0.0317 )
CNPre, CPPre
Real number parameters for α in Eq. (3.18) (Lombardi’s mobility model)
5
for electrons and holes, respectively. (Defaults: CNPre, 1.74 ×10 ,
5
CPPre, 8.84 ×10 )
DELTAN, DELTAP
Real number parameters for δ in Eq. (3.18) (Lombardi’s mobility model)
14
for electrons and holes, respectively. (Defaults: DELTAN, 5.82 ×10 ,
14
DELTAP, 2.05 ×10 )
DIst.epe
Real number parameter for the characteristic length L in computing the
structure-dependent transverse field (Eq. (3.21)), in units of µ .
(Default: 0.5)
EC.Crit
Real number parameter for E crit in Eq. (3.15) (Intel’s local field mobility
4 4
model), in units of V/cm. (Default: 4.2 ×10 for electrons and 3.0 ×10 for
holes)
EC.Unv
Real number parameter for Euniv in Eq. (3.16) (Intel’s local field mobility
5 5
model), in units of V/cm. (Defaults: 5.71 ×10 for electrons and 2.57 ×10
for holes)
ETa
Real number coefficient for η in Eq. (3.22) of Watt’s surface mobility
model. (Defaults: 0.5 for electrons and 0.33 for holes)
EXP.Conc
Real number parameter of c in Eq. (3.6) for analytical low-field mobility
model. (Defaults: see Table 3.3)
EXP.Eper
Real number parameter of β in Eq. (3.15) for Intel’s local field mobility
model (Section 3.1.2.1). (Defaults: 0.5)
EXP1.tem
Real number parameter of b in Eq. (3.6) for the analytical low-field
mobility model. (Defaults: see Table 3.3)
EXP2.tem
Real number parameter of d in Eq. (3.6) for the analytical low-field
mobility model. (Defaults: see Table 3.3)
G.surface
Real number parameter of the reduction factor for the surface mobility
with respect to the bulk one (gsurf in Eq. (C.3)), unitless. (Default: 1.)
Hole
Logical parameter to indicate that the parameters specified in the present
card apply to holes if their values are different for holes from electrons.
These parameters include: EC.Unv, ETa, UCons, UMAx, and UMIn
(Default: false, meaning that parameters apply to electrons)
Int.off
Real number parameter of d in Eq. (3.5) to specify the vertical offset from
the Si/SiO2 interface, in units of µ. (Default: 0.06)
Ncrit
Real number parameter of N0 in Eq. (3.6) for the analytical low-field
mobility model. (Defaults: see Table 3.3).
Qss.conc
Real number parameter of c in Eq. (3.7) to specify the conversion factor
for converting qss into equivalent impurity concentration for mobility
7
reduction, in units of cm-1. (Default: 1.0 ×10 )
UMAx, UMIn
Real number parameters of µ max and µ min in Eq. (3.6) for analytical low-
field mobility model. (Defaults: see Table 3.3)
UConst
Real number parameter of the low-field mobility for general
semiconductor unknown to the program, in units of cm2 V -1 s-1. (No
default)
U0.unv
Real number parameter of µ 0 in Eq. (3.8) for Intel’s surface mobility
model, in units of cm2 V -1 s-1. (Defaults: 783.5 for electrons and 247.4 for
holes)
Vsat
Real number parameter for carrier saturation velocity. (Defaults: see
Table 3.9)
EXAMPLES
4
4.2 ×10 ) in Eq. (3.15) and characteristic length 0.04 µ (default value: 0.03)
in Eq. (3.5).
COMMAND MODELS
The MODELS card sets the temperature for the simulation and specifies
model flags to indicate the inclusion of various physical mechanisms and
models.
SYNTAX
numerical parameters:
PARAMETERS
ABs.coef
The optical absorption coefficient ( α in Eq. (3.43), Section 3.3.2) in units
of cm-1, used in photo-generation (flag PHotogen) modeling.
(Default: 0.0)
ACc.sf
The low-field mobility reduction factor for surface accumulation layers,
used in conjunction with the non-local, transverse-field mobility model
TFldmob (refer to [13] for details). (Default: 0.87)
ANalytic
A logical flag specifying an analytical doping and lattice-temperature
dependent mobility model for silicon only (see Section 3.1.1.3).
(Default: false)
ARora
A logical flag specifying an alternative doping and lattice-temperature
dependent mobility model originally developed by Arora [5] for silicon
and now is available for both silicon and GaAs [6] (see Section 3.1.1.4 for
details). (Default: false)
AUger
A logical flag to specify the use of Auger recombination mechanism (Eq.
(3.47)). (Default: false)
B.Electrons, B.Holes
Real number parameters used in the longitudinal field-dependent mobility
reduction for silicon ( β in Eq. (3.30)). (Defaults: B.Electrons, 1.395;
B.Holes, 1.215)
BGn
A logical flag for including band-gap narrowing mechanism.
(Default: false)
BOltzmann
A logical flag for using Boltzmann carrier statistics. (Default: true)
CCsmob
A logical flag to invoke the mobility model of Dorkel and Leturcq for
silicon [7], which includes carrier-carrier scattering effect, doping and
lattice-temperature dependence (see Section 3.1.1.6). (Default: false)
CONMob
A logical flag to specify use of doping-dependent mobility. The default
model is a table look-up for data in both silicon and GaAs at 300K. When
used with other flags for doping dependent mobility such as ANalytic and
CCsmob, the default model is overritten. (Default: false)
CONSrh
A logical flag to specify use of Shockley-Read-Hall recombination with
doping-dependent carrier lifetimes (Eq. (3.45)). (Default: false).
ENgy.mob
E0
A real number parameter used in the (longitudinal) field-dependent
mobility model for GaAs ( E 0 in Eq. (3.31), Section 3.1.3.2).
3
(Default: 4 ×10 V/cm)
FErmidirac
A logical flag for using Fermi-Dirac carrier statistics. (Default: false)
FLDmob
A logical flag for considering the longitudinal field-dependent mobility
reduction. The default model is Caughey-Thomas formulation [14] (Eq.
(3.30)) for silicon and Thim’s model (with negative differential mobility)
[15] (Eq. (3.31)) for GaAs. (Default: false)
FLUx
The incident photon flux at the y=0 surface (Fph in Eq. (3.43),
Section 3.3.2) in units of photons/cm2. (Real number, default: 0.0)
FMob.new
A logical flag for use of a longitudinal field-dependent mobility model
which is derived from the carrier-temperature dependent model (flag
Hypertang
A logical flag to use hyperbolic tangent form for electron mobility in
GaAs (Eq. (3.32)). This model avoids the negative differential mobility as
exhibited in Thim’s model [15], thus improving the numerical stability.
(Default: false)
IMPAct
A logical flag to include the impact ionization as part of the carrier
generation term in the solution process. The default impact ionization
model is the local field-dependent one (Eqs. (3.37)-(3.38)), and the model
parameters can be changed through IMPACT card. This default model can
be overwritten by one of the following flags: IMP.Engy (IMP.NT),
IMP.NVt, and IMP.JT. (Default: false)
IMP.Engy
A logical flag to indicate that the carrier energy (i.e. temperature)
dependent impact ionization model is to be used as opposed to the local
field dependent model (default of flag IMPAct). There are several carrier-
energy dependent models available (see Section 3.3.1.2 for details) and the
default one is invoked by flag IMP.NT. (Default: false)
IMP.Jt
A logical flag for using carrier-temperature dependent ionization rate (Eq.
(3.39)) in the impact ionization model as described in Eqs. (3.37)-(3.38).
(Default: false)
IMPJt.ratio
IMP.NT
A logical flag for using carrier-temperature dependent ionization model as
described in Eq. (3.42). The unique feature of this model is that the carrier
drift velocity (or equivalently, the current density) does not play any role
in the generation rate due to the impact ionization. (Default: false)
IMP.NVt
A logical flag for using carrier-temperature dependent impact ionization
model in which the carrier (drift) velocity is approximated by the
saturation velocity (Eq. (3.40)). (Default: false)
IMP.TN
A logical flag to indicate that carrier-temperature dependent impact
ionization model is to apply to electrons, whereas the impact ionization
model for holes can be either local-field (no flag IMP.TP) or temperature
(flag IMP.TP set) dependent. (Default: false)
IMP.TP
A logical flag to indicate that carrier-temperature dependent impact
ionization model is to apply to holes, whereas the impact ionization model
for electrons can be either local-field (no flag IMP.TN) or temperature
(flag IMP.TN set) dependent. (Default: false)
INComplete
A logical flag to Indicate that incomplete-ionization of impurities should
be accounted for (Eq. (C.2)). (Default: false)
INV.sf
The low-field mobility reduction factor for surface inversion layer, used in
conjunction with the non-local, transverse-field mobility model TFldmob
(refer to [13] for details). (Default: 0.75)
INTELMOB, INTELMOB.par
Logical flags for using Intel’s local field dependent mobility models (Eqs.
(3.15)-(3.16)). Mobility reduction due to both transverse and longitudinal
fields are taken into consideration simultaneously. The difference between
these two models is in the expression for transverse field reduction.
(Default: false)
Lombardi
A logical flag to specify use of Lombardi local transverse filed dependent
mobility model (Eq. (3.17), also see [8]). (Default: false)
OLdtfld
A logical flag to specify use of the earlier version of UT Austin’s non-local
transverse field dependent mobility model (Eq. (3.24), see also [11]).
(Default: false)
Real number parameters to define the location of the gate oxide region
(left, right, and bottom edges) in conjunction with the use of UT Austin’s
non-local mobility models (flags OLdfld and TFldmob). (No default
values but initialized to 0.0)
PHotogen
PRint
A logical flag to print the status of all models and a variety of coefficients
and constants. (Default: false)
SRFmob
A logical flag to invoke Watt’s surface mobility model (Eqs. (3.22)-
(3.23)). (Default: false)
SRH
A logical flag to specify the Shockley-Read-Hall recombination (Eq.
(3.44)) with constant carrier lifetimes. (Default: true)
STrfld
A logical flag to specify that the surface transverse field is evaluated by the
structure information (Eq. (3.21)) rather by definition (Eq. (3.20). This
flag is usually used together with one of the transverse field dependent
mobility models such as INTELMOB. (Default: false)
TAU.WN, TAU.WP
TEmperature
Specifies the environment temperature, in units of Kelvin. (Default: 300)
TFldmob
A logical flag to invoke the second version (as opposed to flag OLdfld) of
UT Austin’s non-local transverse field mobility model (Eq. (3.28), see also
[13]). It is based on an extended version of the Schwarz-Russek
formulation [12]. (Default: false)
User1
A logical flag to specify Arora’s (ARora) doping and lattice temperature
dependent mobility model but with the user-definable parameters (refer to
Section 3.1.1.5). (Default: false)
EXAMPLES
ionization model IMP.JT for the electron system is specified. The energy
relaxation time ratio in computing impact ionization rate is set to be 0.81.
COMMAND OPTIONS
SYNTAX
PARAMETERS
CPUFile
Character string for the name of the file to log the CPU statistics during
PISCES run. (Default: null)
CPUStat
A logical flag to indicate that CPU statistics during PISCES run is to be
outputted to the file specified by CPUFile. (Default: false)
If PSfile is used, a PostScript file named print.ps will be generated for the
plot and can later be printed to obtain a hard copy using laser printers.
Plots will be scaled to the size of the specific device. Also note that on
color graphics terminals, the different line types are implemented as
different colors; on the black and white monitors they are implemented as
dot and line patterns. (Default: login terminal)
PLOTFile
A character string for name of the output file which is generally defined by
the plot device. For example, a graphics terminal will use the terminal as
the output file. Printers may have the output file be a spooler. The graphics
output file can explicitly be set by the PLOTFile command. All graphics
output will then be routed to the given file. Note that the contents of the
X.Offset, Y.Offset
Real number parameters, in units of inches, to indicate the x- and y-offsets
from the bottom-left corner of the screen, respectively. (Defaults: 0.0)
X.Screen, Y.Screen
Real number parameters, in units of inches, for the physical width and
height of the screen, respectively. They are set automatically, depending
on the plot device. They, however, can be altered for special effects (split
screen plots, for instance). (Defaults: device’s size)
EXAMPLES
COMMAND PLOT.1D
The PLOT.1D card plots a specific quantity along a line segment through
the device (mode A), or plots an I-V curve of data (mode B). In mode A,
one of the solution variables is plotted versus distance into the device. For
vector quantities, the magnitude is plotted. In mode B, terminal
characteristics can be plotted against each other by choosing the value to
be plotted on each axis.
SYNTAX
PARAMETERS
ABsolute, X.ABsolute
Logical flags to indicate use of absolute values for quantities to be plotted
in y- and x-axes, respectively. (Defaults: false)
AScii
Specifies that the output file will have an ascii format (the default binary
format is for use with the Stanford dplot system). (Default: false)
BAND.Val
A logical flag for plotting valence-band potential. (Default: false)
BAND.Con
A logical flag for plotting conduction-band potential. (Default: false)
DAta
Logical flag to indicate an HP measured I-V file will be used as an input
file for plotting. (Default: false)
DOping
A logical flag for plotting doping. (Default: false)
E.field
A logical flag for plotting electric field. (Default: false)
ELectrons
A logical flag for plotting electron concentration. (Default: false)
GEN.Elec, GEN.Hole
Logical flags to plot the generation rates due to impact ionization caused
by electron and hole currents, respectively. (Defaults: false)
Holes
A logical flag for plotting hole concentration. (Default: false)
IMpact (II.gener)
A logical flag for generation rate due to impact ionization (G in
Eq. (3.37)). (Default: false)
INFile
Character string for name of the input file in which the solution is stored.
(Default: null)
INTegral
Plots the integral of the specified ordinate. (Default: false)
LIne.type
Specifies the line type for the plotted curve. (Default: 1)
MIn.value, MAx.value
Specify minimum and maximum values, respectively, for the ordinate of
the graph. Default: found automatically from the data to be plotted.
NEGative
Negates the ordinate values. PISCES-II by default will order the plot
coordinates by abscissa value; this ordering will result in unusual plots for
IV curves with negative resistance, for example. (Default: false)
NET.CArrier
Plot net carrier concentration. (Default: false)
NET.CHarge
Plot net charge concentration. (Default: false)
NO.Axis
Indicates that the axes for the graph are not to be plotted. (Default: false)
NO.Clear
Indicates that the screen is not to be cleared before the current plot so that
several curves can be plotted on the same axis. (Default: false)
NO.Order
A logical flag to force PISCES to plot the data points as they naturally
appear. (Default: false)
Outfile
Character string for name of an output file in which the data points plotted
are to be put. If specified, the graphics output will be directed to that file.
For further discussion, see the Options card. (Default: from Options
card)
PAuse
Causes PISCES-II to stop at the end of the plot so that a hard copy may
be made before continuing. Execution can be resumed by hitting a
carriage return. (Default: false)
POInts
Marks the data points on the plotted curve. (Default: false)
POTential
Plot mid-gap potential. (Default: false)
QFN, QFP
Logical flags to plot electron and hole quasi-fermi levels, respectively.
(Defaults: false)
Recomb
Plot net recombination. (Default: false)
Spline, NSpline
The Spline option indicates that spline-smoothing should be performed on
the data using NSpline interpolated points (maximum is 500). (Default:
Spline, false; NSpline, 100)
Unchanged
A synonym for NO.Axis and NO.Clear, but additionally it forces the use
of the previous axis bounds so that a number of curves can easily be put on
the same axis. (Default: false)
VAcuum
Logical flag to plot the vacuum energy level, a useful quantity for
simulation of heterostructures because it is always continuous even at the
abrupt material interface. (Default: false)
VELO.Ele, VELO.Hol
Logical flags to plot the electron and hole velocities, respectively.
(Defaults: false)
X.AXis, Y.Axis
Character strings to indicate which quantities are to be used as x- and y-
axes, respectively.
X.Component, Y.Component
Logical flags to force the x or y components, respectively, of any vector
quantities to be plotted, as opposed to the default magnitude. (Defaults:
false)
For rapidly varying quantities, the use of logarithmic values is often more
revealing. Since many of the quantities may become negative, PISCES
actually uses the following expression
X.MAx, X.MIn
Real number parameters to specify the maximum and minimum values for
the abscissa to be plotted. (Defaults: the maximum and minimum abscissa
values in the data to be plotted)
XMole
A logical flag to plot the mole fraction of a compound material. For
quaternaries, only x nor y mole fraction is plotted. (Default: false)
EXAMPLES
COMMAND PLOT.2D
SYNTAX
PARAMETERS
Crosses
A logical flag to plot crosses at the locations of grid points. (Default: false)
Boundary
A logical flag to indicate that the boundaries around the device and
between regions are to be plotted. (Default: false)
Depl.edge
A logical flag to indicate that depletion edges are to be plotted. Note that
depletion edges can only be plotted when a solution is present. (Default:
false)
Grid (Mesh)
A logical flag to plot the grid including lines which delineate elements.
(Default: false)
Flip.x
A logical flag to flip the plot about the y-axis; i.e., it negates all x
coordinates so that the plot is mirrored. (Default: false)
Junction
A logical flag to specify that the junctions from the doping profile are to be
plotted. (Default: false)
LAbels
A logical flag to indicate that room is to be made for color contour labels
on the right side of the plot device. (Default: false)
NO.Clear
A logical flag to specify that the screen is not to be cleared before plotting.
(Default: false)
NO.Diag
A logical flag to indicate not to plot the diagonals for a rectangle-based
mesh. (Default: false)
NO.Fill
A logical flag to force the program to draw the device area plotted to scale.
If this option is not specified, the plot will fill the screen and the triangles
will appear distorted. (Default: false)
NO.TIc
A logical flag to indicate that tic marks are not to be included around the
plotted area. (Default: false)
NO.TOp
A logical flag to indicate that tic marks are not to be put on the top of the
plotted region. (Default: false)
Outfile
Character string for name of an output file in which the data points plotted
are to be put. If specified, the graphics output will be directed to that file.
For further discussion, see the Options card. (Default: from Options
card)
Pause
A logical flag causing the program to stop at the end of the plotting so that
a hard copy may be made before continuing. Execution can be resumed by
hitting a carriage return. (Default: false)
EXAMPLES
COMMAND PRINT
The PRINT card prints specific quantities at grids within a defined area of
the device.
SYNTAX
PARAMETERS
Current
A logical flag to print currents (electron, hole, conduction, displacement,
and total) at each grid for the present solution. (Default: false)
Elements
A logical flag to print information on triangular elements (sequence
number, nodes, and material). (Default: false)
Geometry
A logical flag to print geometrical information (vertex coordinates) of
triangles. (Default: false)
Material
A logical flag to print material information (permittivity, band-gap, etc.)
including the value of the doping dependent mobility and lifetime (if
model flags specified) at each grid. (Default: false)
P.CURR1, P.CURR2
Logical flags to print currents (electron, hole, conduction, displacement,
and total) at each grid for previous first and second solutions, respectively.
(Defaults: false)
POints
A logical flag to print grid information (coordinates, doping, etc.).
(Default: false)
P.QUE1, P.QUE2
Logical flags to print space charge, recombination rate, and electric field
for the previous first and second solutions, respectively. (Defaults: false)
P.SOL1, P.SOL2
Logical flags to print previous first and second solutions (ψ, n, p, and
quasi-Fermi potentials), respectively. (Defaults: false)
Que
A logical flag to print space charge, recombination rate, and electric field
for the present solution. (Default: false)
Solution
A logical flag to print the present solution (ψ, n, p, and quasi-Fermi
potentials). (Default: false)
X.Compon, Y.Compon
Logical flags to specify how any of the various vector quantities (current
and field) should be printed. The default is the magnitude of the vector.
X.Compon specifies that the magnitude of the x-component of the vector
be printed, while Y.Compon specifies the y-component. Only one of
these (or neither) can be specified on a single card. (Defaults: false)
EXAMPLES
COMMAND REGION
The REGION card assigns material type and its composition, i.e., mole
fraction profile for a compound semiconductor, to a rectangular region
bounded by high and low node numbers in x- and y-directions. Every
triangular element must be given a certain type of material. If there is an
overlap among region assignments, the latter one takes precedence.
Currently, there are four types of compound semiconductors available in
the code and they are GexSi1-x, AlxGa1-xAs, AlxIn1-xAs, and GaxIn1-
xAsyP1-y. The composition profile can be either constant or linear.
SYNTAX
NUmber = <integer>
position:
mole spec:
PARAMETERS
ALGaas
A logical flag indicating an AlxGa1-xAs region. (Default: false)
ALInas
A logical flag indicating an AlxIn1-xAs. (Default: false)
GAAs
A logical flag to indicate a gallium arsenide (GaAs) region.
(Default: false)
GAInasp
A logical flag indicating a GaxIn1-xAsyP1-y region. (Default: false)
INsulator
A logical flag to indicate a general insulator region. (Default: false)
NItride or SI3n4
A logical flag to indicate a nitride (Si3N4) insulator region. (Default: false)
NUmber
This parameter selects the rectangular region in question. Currently in the
code the maximum number of regions is set to 1,000, but can be changed
under the user's request.
Oxide or SIO2
A logical flag to indicate a silicon dioxide (SiO2) region. (Default: false)
SApphire
A logical flag to indicate a sapphire (Al2O3) region. (Default: false)
SEmiconductor
A logical flag to indicate a general semiconductor region. Forty material
parameters can be defined by the user through MATERIAL card. (Default:
false)
SIGe or GEsi
A logical flag indicating a GexSi1-x region. (Default: false)
SILicon
A logical flag to indicate a silicon (Si) region. (Default: true)
STrained
A logical flag to indicate that the material is a strained layer. Currently in
the code, it is only applied to GaxIn1-xAsyP1-y and only the bandgap and
the effective mass for holes (both light and heavy ones) are affected by this
flag (for details see Section 4.4.3). (Default: false)
XEnd – XIinital
mole ( d ) = XIinital + ----------------------------------------- ⋅ ( d – d low ) (M.4)
d high – d low
XMole, YMole
These real number parameters are used to specify the mole fraction for
compound semiconductors, be binary (GexSi1-x), ternary (AlxGa1-xAs and
AlxIn1-xAs), or quaternary (GaxIn1-xAsyP1-y). If the composition profile is
flat (constant) in the region, the mole fraction is specified using XMole for
x in binary/ternary and XMole and YMole for x and y, respectively, in
quaternary.
YLinear
A logical flag to indicate that the linear composition profile changes along
the y-direction instead of x-direction. (Default: false)
EXAMPLES
COMMAND REGRID
SYNTAX
files:
PARAMETERS
ABsolute
A logical flag to specify that the absolute value of the quantity is to be
used. (Default: false)
AScii
A logical flag to indicate that all mesh files and triangle trees (not
including DOPFile) for this card should be done in ASCII rather than
binary (a default). (Default: false)
CHange
A logical flag to determine whether to use the magnitude (false) or the
difference (true) of a variable in a triangle as the criterion of refinement.
(Default: true)
COs.angle
A real number for the cosine of the angle, which defines the “obtuse
criterion'' to limit the creation of obtuse angles in the mesh. If regrid would
create a triangle with an angle whose cosine ( cos ) is less than -
COs.angle, nodes are added so that this does not occur. The test can be
turned off locally by using the Ignore flag in this card. It can be turned off
everywhere by using a value of COs.angle greater than 1.0. The default
is to turn it off everywhere. (Default: 2.0)
DOPFile
Character string for the name of the binary mesh file, which contains the
doping for the device (see DOPING card). If set to a non-null string,
interpolation of doping values at any newly created grid points to re-dope
the structure based on the initial doping specification (a default action)
will not occur. (Default: null)
DOPIng
A logical flag to indicate regriding based on doping concentration (in units
of cm-3). (Default: false)
ELEctrons
A logical flag to indicate regriding based on electron concentration (in
units of cm-3). (Default: false)
EL.field
A logical flag to indicate regriding based on electric filed (in units of V/
cm). (Default: false)
Holes
A logical flag to indicate regriding based on hole concentration (in units of
cm-3). (Default: false)
Ignore
An integer number for region(s). In these regions, regriding will not
happen (i.e., ignored), nor will be smoothed after regriding. This flag is
similar to REGion in function, but opposite in effect. (Default: 0, i.e.,
none ignored)
LOCaldop
A logical flag related to the minority carrier regrids. If set, when the
minority carrier concentration exceeds the local doping, the grid will be
refined. (Default: false)
LOGarithm
A logical flag to use logarithmic instead of the actual value of the quantity
concerned. If set, parameter STep (RAtio) is interpreted in the logarithm.
(Default: false.)
MAx.level
An integer for the maximum level of any triangle relative to the original
mesh. It defaults to one more than the maximum level of the grid, but can
be set to a smaller value to limit refinement. Values less than or equal to
zero are interpreted relative to the current maximum level. (Default:
dynamic)
MIn.carr
A logical flag to indicate regriding based on minority carrier concentration
(in units of cm-3). (Default: false)
NET.CHarge
A logical flag to indicate regriding based on the net charge density (in
units of cm-3). (Default: false)
NET.CArr
A logical flag to indicate regriding based on net carrier concentration (in
units of cm-3). (Default: false)
Outfile
Character string for the binary output mesh file, which is necessary if the
mesh is to be used for subsequent runs. A history of the triangle tree is
always generated to assist further regriding steps and its name is the one
specified by Outfile and concatenated by the letters “tt” to the end. In the
case that this file is used as an input mesh file for a PISCES run including
a regrid action, the program will look for a file with the same name as the
input mesh file plus “tt” at the end. (Default: null)
Potential
A logical flag to indicate regriding based on electrostatics potential which
usually refers to the intrinsic Fermi level (in units of V). (Default: false)
QFN, QFP
Logical flags to indicate regriding based on electron and hole quasi-Fermi
levels, (in units of cm-3), respectively. (Defaults: false)
Region
An integer for region(s) in which mesh is to be refined according to the
user criterion. Note that other regions might be refined as well as a side
effect to maintain well-shaped triangles. (Default: all regions)
SMooth.key
An integer parameter which has the same meaning as SMooth.key in
MEsh card. See MEsh card for details. (Default: 0)
STep (RAtio)
Real number parameter as the numerical criterion for refining a triangle. It
is interpreted in logarithm if flag LOGarithm is set. (No default)
EXAMPLES
COMMAND SOLVE
The SOLVE card instructs PISCES to perform a solution for one or more
specified bias points.
SYNTAX
ac:
PARAMETERS
AC.analysis
A logical flag to indicate that AC sinusoidal small-signal analysis be
performed after the DC condition is solved for. Note that the full Newton
method (2 carriers) must be used for this analysis. Use of flag G.Debug
in the Options card will provide some detailed information on the AC
solution procedure. (Default: false)
AScii
A logical flag to indicate that the solution file saved in Outfile will be
ASCII as opposed to binary (a default). (Default: false)
BAnd
A logical flag to include the band profile (vacuum level, conduction and
valence band edges) in the output solution file in addition to potential,
carrier concentrations, etc. (Default: false)
Currents
A logical flag to indicate that the electron, hole, and displacement
currents, and the electric field, will be computed and stored with the
solution. (Default: false)
Dt or TSTEp
Real number for the time-step to be taken during time transient analysis.
For automatic time-step runs (flag TAuto in MEthod card), it is used to
select the first time step only. (Default: 0.0)
ELectrode
An integer for the number of the electrode being stepped. If more than one
electrode is to be stepped, ELectrode should then be an n-digit integer,
where each of the n-digits is a separate electrode number. Note that if there
are 10 electrodes, don't put electrode 0 first in the sequence!). (Default: 0)
EXtrapolate
A logical flag to indicate that the initial guess for the solution at the
present bias is to be obtained by linear extrapolation of the previous two
solutions which are either just obtained in the current PISCES run or
loaded in using parameters IN1file and IN2file in LOad card.
(Default: false)
FRequency
Real number parameter for the frequency, in units of Hz, at which AC
analysis is to be performed. (No default)
INitial
A logical flag to indicate the solution at the thermal equilibrium (i.e., zero
bias) is to be sought. Note that the first bias point in any PISCES run must
be set with this flag unless there is a solution loaded using LOad card.
(Default: false)
IStep
Real number parameter for the current increment to be added to one or
more electrodes, as specified by the integer assigned using parameter
ELectrode in this card. (Default: 0.0)
LOCal
A logical flag to use the carrier quasi-Fermi levels in the previous solution
(INFile in LOAd card) as the initial guess. (Default: false)
LOWfreq.ac
A logical flag to use zero-frequency for AC analysis (see Section 5.1.4 for
details). (Default: false)
MAx.inner
Integer parameter for the maximum number of SOR iterations.
(Default: 25)
N.bias, P.bias
Real number parameters to set the levels for electron and hole quasi-Fermi
potentials, respectively. These parameters only apply to those carriers
which are not being solved for in the PISCES run (e.g. holes when only
electrons are specified to solve in Solve card). If N.bias or P.bias is not
specified, then program will either choose local quasi-Fermi potentials
based on bias and doping (see Chapter 2 of [3]) or if Fix.qf is set on the
MEthod card, set the quasi-Fermi levels where applicable to values which
produce the least amount of free carriers (maximum bias for electrons and
minimum bias for holes).
NSteps
The number of bias increments (steps) to be taken; i.e., if VStep (IStep)
is specified, the specified electrode is incrementted NSteps times.
(Default: 0.0)
Outfile
Character string for the name of the binary output file to save the solution
at this bias point. If an electrode is stepped so that more than one solution
is generated by this card, the last non-blank character of the supplied file
name will have its ASCII code incrementted by one for each bias point in
succession, resulting in a unique file for each bias point. (Default: null)
PREvious
A logical flag to specify that the solution at the previous bias point, either
loaded in using LOad card or just solved, be used as the initial guess for
the solution at the present bias point. (Default: true)
PROject
A logical flag to indicate Newton Projection method (NPM) is to be used
in obtaining the initial guess to the solution at the present bias based on the
previous solution. For details of NPM, see Section 5.4. (Default: false)
Qmr
A logical flag to use QMR (Quasi-Minimal Residual) method for AC
analysis (see Section 5.1.3 for details) (Default: true)
Ramptime, ENdramp
Real number parameters to specify the time duration of a bias ramp and
the ending instant of the ramp, respectively, for transient analysis, in units
of seconds. Note that bias ramp must be linear. Specifically, is the ramp
begins at t = t0, then it ends either at t = t0 + Ramptime or at t =
ENdramp, depending on which parameter is used. (Defaults: 0.0)
S.omega
A real number parameter for the relaxation factor used in SOR method ( ν
in Eqs. (5.13)-(5.14)). Note that it is not a frequency. (Default: 1.0)
SOr
A logical flag to use SOR method for AC analysis (refer to Section 5.1.1.
and [36] for details) (Default: false)
TErminal
An integer parameter for the sequence number of contact(s) to which the
AC bias (i.e., stimulus) is to be applied. More than one contact number
may be specified (via concatenation), but each will be solved separately.
Each contact that is specified yields a column of the admittance matrix.
(Default: all contacts)
TOlerance
Real number parameter used as the criterion for the iteration to terminate
–5
in AC analysis, unitless. (Default: 1 ×10 )
TSTOp (TFinal)
Real number parameter to specify the end of the time interval to be
simulated, in units of seconds. Thus the simulation begins at t = t0, it will
end at t = TSTOp (TFinal). Alternatively, NSteps can be used to signal
the end of the interval; That is, the final time would be t = t0 + NSteps x
TSTEp. (Default: 0.0)
VSS
Real number parameter for the magnitude of the applied small-signal bias
(stimulus) (Vi in Eq. (2.19) of [4]). (Default: 0.1 x kT/q where T is the
environment temperature)
VSTep
Real number parameter for the voltage increment to be added to one or
more electrodes, as specified by the integer assigned to parameter
Electrode in Solve card. (Default: 0.0)
EXAMPLES
Bias point # V1 V2 V3
1 0.0 0.5 0.5
2 1.0 0.5 0.0
3 2.0 0.5 0.0
4 3.0 0.5 0.0
5 4.0 0.5 0.0
6 5.0 0.5 0.0
The solutions for these bias points will be saved to the files OUT1, OUTA,
OUTB, OUTC, OUTD and OUTE. Note that the initial guess for the first
bias point is obtained directly from the preceding solution because the
PREvious option was specified. The initial guesses for bias points 2 and
3 will also be obtained as if PREvious had been specified since two
electrodes (numbers 1 and 3) had their biases changed on bias point 2.
However, for bias points 4, 5 and 6, PISCES-II will use a projection to
obtain an initial guess since starting with bias point 4, both of its preceding
solutions (bias points 2 and 3) only had the same electrode bias (number 1)
altered.
ramp is at t = 45ns). The device is then solved at this bias for another 55ns
(out to 100ns). Note that again each solution is saved in a separate file
(DOWN1, DOWN2, etc.) and that no initial time-step was required since
one had been estimated from the last transient solution for the previous
SOLVE card. Finally, the fourth SOLVE card performs the steady-state
solution at V1 = -1 and V2 = 0.0
COMMAND SPREAD
SYNTAX
Width = <real>
Upper = <logical> LOwer = <logical>
specifics:
PARAMETERS
Encroach
A real number factor which defines the abruptness of the transition
between a distorted and non-distorted mesh. The transition region
becomes more abrupt with smaller Encroach factors (the minimum is
0.1). An important note: depending on the characteristics of the
undistorted mesh, very bad triangles (long, thin, and obtuse) may result if
Encroach is set too low. (Default: 1.0)
GRAding
A real number parameter to specify a grid ratio (identical to the Ratio
parameter on the X.Mesh and Y.Mesh cards) to produce a non-uniform
grid in the distorted region. (Default: 1.0)
LEft, Right
Logical flags to specify that the left and right sides of the grid,
respectively, will be distorted. (Defaults: false)
Upper, LOwer
Integer parameters to specify the upper and lower y-grid lines,
respectively, between which the distortion will take place. (No default)
Vol.ratio
Real number parameter as the ratio of the downward displacement of the
lower grid line to the net increase in thickness. It is ignored if Y.Lower is
specified. (Default: 0.44 which is derived from oxide/silicon growth/
consumption ratio during oxidation such that the Si/SiO2 interface can
move correctly)
Width
A real number parameter to specify the width from the left or right edge
(depending on the LEft and Right parameters) of the distorted area. The
actual x-coordinate specified by Width (min[x] + Width for LEft and
max[x] - Width for Right) will lie in the middle of the transition region
between the distorted and undistorted grid regions, in units of µ . (No
default)
EXAMPLES
COMMAND SYMBOLIC
SYNTAX
Carriers = <integer>
{ Electrons = <logical> | Holes = <logical> }
ENGY.Elect = <logical> | ENGY.Hole = <logical>
ENGY.Lat = <logical>
options:
PARAMETERS
Carriers
ELectrons, Holes
Logical flags to indicate which type of carriers is to be solved if Carriers
is set to 1. (Defaults: false)
ENGY.Elect, ENGY.Hole
Logical flags to indicate if electron or hole or both carrier energy balance
equation(s) is to be solved in the simulation. (Defaults: false)
ENGY.Lat
Infile, Outfile
Character strings (up to 20 characters long) to specify the input and output
(binary) file names, respectively, for the symbolic factorization. Note that
these binary files can be quite large, so it is advised not to use this feature.
(In some computing environments it is also faster to compute the symbolic
information for each run than to read it from the file). (No default, strings
set to null)
Min.degree
A logical flag to specify the use of a minimum degree ordering of the
pivots for decomposition in order to reduce the size of the generated L and
U matrices, hence to reduce the amount of CPU time spent in solving
linear systems. This parameter is definitely recommended. (Default: true)
Newton, Gummel
Logical flags to specify the Newton simultaneous or Gummel sequential
method in solving the Poisson’s, electron and hole continuity equations. In
the case of solving the Poisson’s equation only (Carriers = 0), Gummel
method is the default one. When one of or both carrier energy balance
equations is to be solved, Newton method must be used for the current
implementation. (Defaults: false)
Print
A logical flag to indicate that information about the memory allocated for
the run is to be printed to the standard output file. (Default: false)
Strip
A logical flag to specify redundancy (zero couplings) be removed from the
symbolic map, and is naturally on. (Default: true)
EXAMPLES
COMMAND TITLE
SYNTAX
Title
<character string>
PARAMETERS
None.
EXAMPLES
COMMAND VECTOR
The VECTOR card plots vector quantities over an area of the device
defined by the previous PLOT.2D card.
SYNTAX
PARAMETERS
Clipfact
A threshold below which vectors are not plotted. (Default: 0.1)
E.field
Plot the electric field. (Default: false)
LIne.type
An integer parameter for the vector line type for plotting. (Default: 1, a
solid line)
LOgarithm
A logical flag to specify use of logarithmic value of the magnitude. By
default, all vectors in a linear plot are scaled by the maximum magnitude
of the quantity of interest over the grid, while if LOgarithm is set, the
default scaling uses the minimum (non-zero) magnitude. (Default: false)
MAximum, MInimum
Real number parameters for the minimum and maximum values of the
quantities to be plotted. These values will be printed during the execution
of a plot. With these parameters is possible to plot two bias conditions or
two devices with the same scaling. (Defaults: 0.0)
Scale
A real number parameter used as the scale factor by which all magnitudes
are multiplied. (Default: 1.0)
EXAMPLES
The X.MESH and Y.MESH cards specify the location of lines of nodes
in a rectangular mesh.
SYNTAX
Node = <integer>
location:
Location = <real>
ratio:
Ratio = <real>
PARAMETERS
Location
Real number parameter to specify where to locate the line, in units of µ .
(No default)
Node
An integer parameter for the numbering of the line in the mesh. Lines
should be assigned consecutively, beginning with the first and ending with
the last. (No default)
Ratio
Real number parameter as the ratio used for the interpolation of lines
between the given lines, unitless. The spacing grows/shrinks by ratio in
each subinterval, and the ratio should usually lie between 0.667 and 1.5.
(Default: 1.0)
EXAMPLES
Curve Tracer
This project was supported by Semiconductor Research Corporation Contract Numbers, SRC 94-SJ-
116 and SRC 93-MC-106. We would also like to acknowledge the mentorship of Dr. Ben Tseng of
AMD as well as AMD’s support of both a summer internship and ongoing access to their facilities and
device technology during the period of this study.
In Technology CAD, the use of software to simulate the testing of semiconductor devices is known as
the Virtual Instrumentation. A virtual instrument should be able to automatically generate simulation
data, e.g., I-V curve along a bias sweep, given only the simple specifications a user would input to a
real programmable instrument testing an IC device. Numerical device simulators such as PISCES
provide a means of creating virtual devices and simulating electrical tests on the devices. However,
these simulators cannot trace through I-V curves with sharp turns unless the user carefully controls the
bias conditions near these turns -- a tedious and time-consuming process. This deficiency prompted
the creation of Tracer.
Tracer is a C program which automatically guides PISCES and other semiconductor device
simulators through complex I-V traces and is ideally suited for device-failure phenomena such as
latchup, BVCEO, and electrostatic-discharge (ESD) protection. Given a PISCES input deck and a
specification file with a PISCES-like syntax, a simulation can be run over any current or voltage range
without the user’s intervention. Tracer is limited to DC, one-dimensional traces, i.e., only one
electrode can be swept per run. It sweeps the electrode by dynamically setting the most stable bias
condition at each solution point. Additionally, Tracer has the ability to maintain zero-current bias
conditions at one or two electrodes during the trace, even at low device-current levels where such bias
conditions are unstable using traditional device simulation. The theory implemented in the Tracer
program can be found in [1].
where
• inputfile is the name of the PISCES input deck which defines the device structure to be
simulated and specifies what physical models are to be used. Basically, it contains everything
in a normal PISCES deck except the solve card specifications (SECTION 10).
• outputfile is an optional specification of the name of the file where the simulation data is
to be written (SECTION 11). If outputfile is not given, the name of the output file defaults
to inputfile.out.
The trace specification file, tracefile, is similar to a PISCES or SUPREM input deck. Each line
begins with a word designating what type of statement (or “card”) it is. The four available cards are
CONTROL, FIXED, OPTION, and SOLVE. Also, a line may start with a “$” for comments. Such
lines are ignored during program execution. The cards may appear in any order. A card may be
continued on following lines by placing a “+” at the beginning of each subsequent line. The “+” should
be separated from the parameters on the line by at least one space.
Each option in a card should have the following structure: “param = paramvalue”. Spaces
separating the “=” sign are optional. The parameters for each card are described below. As with
PISCES syntax, parameter names and values are not case-sensitive and may be abbreviated provided
they remain unambiguous. Square brackets, [], enclosing a parameter indicate that it is optional (note
that some of these parameters are optional only in the sense that they will default to a certain value if
not specified in tracefile). A vertical line, |, represents a logical OR -- only one of a list of
parameters separated by “|” signs can be specified.
All electrodes in the device must have representation in the tracefile. Each electrode
must appear as one, and only one, of the following: the CONTROL electrode, a FIXED electrode, or
an open contact (OPENCONT1 or OPENCONT2) on the SOLVE card.
9.1.2 Syntax
NUM=<int> CONTROL=<char> [BEGIN=<real>] [INITSTEP=<real>]
[ENDVAL=<real> | STEPS=<int>]
9.1.3 Parameters
• NUM is the sequence number of the electrode in the PISCES input deck designated as the
control electrode, whose voltage or current is swept through the trace. Its integer value must be
between 1 and 9, inclusive. (Default: none)
• CONTROL is either VMAX, IMAX, or STEP. VMAX denotes that a maximum voltage on
the control electrode, specified by ENDVAL, is used as the upper bound on the trace. IMAX
denotes that ENDVAL specifies a maximum control-electrode current for the trace. STEP
signifies that the trace will proceed for a certain number of simulation points, specified by the
STEPS parameter. In most cases VMAX or IMAX will be used because it is not known how
many simulation steps will be taken to reach a certain voltage or current. (Default: none)
• BEGIN is the value of the voltage, in volts, at the starting point of the curve trace for the
electrode designated by NUM (the control electrode). If an initial solution is performed by
Tracer, BEGIN should be 0.0. If a previous solution is loaded into the input deck at the start
of Tracer (see SOLVE card below), BEGIN should be equal to the voltage of the control
electrode in this solution. (Default: 0.0V)
• INITSTEP is the initial voltage increment, in volts, of the control electrode. Thus, at the
second solution point the control electrode will have a voltage of BEGIN + INITSTEP. A
recommended initial step size is 0.1V. The sign of INITSTEP determines the direction in
which the curve tracing will initially be proceeded. If INITSTEP proves to be too large and
PISCES cannot converge on the second solution point, Tracer will automatically reduce
INITSTEP until convergence is attained, then proceed with the trace from this point. (Default:
0.1V)
• ENDVAL is used when CONTROL=VMAX or IMAX. Tracer stops tracing when the
voltage (CONTROL=VMAX) or current (CONTROL=IMAX) of the specified electrode
equals or exceeds the value specified by ENDVAL. Note that it is the absolute values of the
voltage or current and of ENDVAL which are compared. (Default: 10.0V
(CONTROL=VMAX), 10.0A/µm (CONTROL=IMAX))
• STEPS is used when CONTROL=STEP. It specifies the number of solution points Tracer
should find. (Default: 10)
9.1.4 Examples
1. Electrode 3 is the control electrode. Tracer will initially proceed in the negative-voltage direction
with an initial step of -0.1V. Tracer will proceed until the absolute value of the control current
equals or exceeds 3A/µm.
2. Electrode 4 is the control electrode. Tracer will run until 65 solutions are found, starting at
v4=0.0V with an initial v4 step of 0.5V.
electrodes in a simulation, in inputfile the user must create contact cards with the “current” option
for each of these electrodes (see SECTION 10).
9.2.2 Syntax
NUM=<int> [TYPE=<char>] [VALUE=<real>] [RECORD=<char>]
9.2.3 Parameters
• NUM is the sequence number of an electrode in the PISCES deck. Its integer value must be
between1 and 9, inclusive. (Default: none)
• TYPE is either VOLTAGE or V for a voltage source or CURRENT or I for a current source.
(Default: VOLTAGE)
• VALUE is the fixed value of the current or voltage for the electrode specified by NUM.
VALUE has units of either volts or amps/µm, depending on the specification of TYPE. Note
that the specification of VALUE is optional since it is merely for reference and is not used by
Tracer. (Default: 0.0)
• RECORD is either YES or NO. For RECORD=YES, the simulated current is recorded in the
output file for a fixed-voltage electrode, while the simulated voltage is recorded for a fixed-
current electrode. (Default: NO)
9.2.4 Examples
In every Tracer solution, electrode 1 has a voltage of 0.0V. The current in this node is recorded in
outputfile at every solution point.
9.3.2 Syntax
Simulations with one or two open contacts:
9.3.3 Parameters
• ABSMAX is the maximum current allowed in an open contact and is only relevant when open
contacts are used and voltage biases are applied to these contacts. Convergence is satisfied
when either the ABSMAX or RELMAX condition is met. (Default: 1.0 ×10 A/ µ m)
– 19
• DAMP is a number between 0 and 1.0 determining how quickly Tracer will converge on an
open-contact solution using voltage biasing. The closer DAMP is to 1.0, the more quickly
Tracer will converge, but there is also an increased chance of slower convergence due to
overshoot. Usually the user should not be concerned with the value of DAMP. (Default: 0.9)
• TRYCBC is used only if there is an open contact. Tracer will only attempt to use zero-
current biasing when the current of the control electrode is greater than TRYCBC. Otherwise,
voltage biasing is used. In most cases the user does not have to worry about this parameter.
– 17
(Default: 1.0 ×10 A/m)
• ANGLE1, ANGLE2, and ANGLE3 are critical angles (in degrees) affecting the smoothness
and step size of the trace. They are described in detail in [1]. If the difference in slopes of the
last two solution points is less than ANGLE1, the step size will be increased for the next
projected solution. If the difference is between ANGLE1 and ANGLE2, the step size remains
the same. If the difference is greater than ANGLE2, the step size is reduced. ANGLE3 is the
maximum difference allowed, unless overridden by the MINDL parameter. ANGLE2 should
always be greater than ANGLE1 and less than ANGLE3. (Defaults: ANGLE1 = 5˚,
ANGLE2 = 10˚, ANGLE3 = 15˚)
• ITLIM is the maximun number of Newton loops for a given solution as specified in the
method card of the PISCES input deck. The user should make sure that the value of ITLIM
specified here is the same as that in the input deck. In certain cases, a PISCES solution may be
aborted in Tracer because the solution will not converge within the given number of
iterations. In some of these cases Tracer will try to redo the solution with a doubled number
of iterations. If ITLIM is specified here, such attempts will be made. If there is no itlim
statement or ITLIM=0, no attempts will be made. It is recommended that ITLIM be set to a
low value, around 10 or 15 (or at least high enough to allow convergence of the initial
solution). However, for GaAs devices a larger ITLIM of 20 or 25 is recommended.
(Default: 0)
• MINCUR is the value of the control current, in A/µm, above which Tracer carefully controls
step size and guarantees a smooth trace. Below this current level, the program simply takes
voltage steps as large as possible, i.e., as long as numerical convergence can be achieved,
without regard for smoothness. If MINCUR is set to 0.0, Tracer will not begin smoothness
control until it is past the first sharp turn in the I-V curve. This value should be used when the
user is only interested in the rough location of a break in the curve, such as the breakdown
voltage of a single-junction device. If smoothness is required, a lower value should be
specified. Setting MINCUR below 1.0 ×10 A/ µ m is not recommended because Tracer has
– 15
• MINDL is the minimum normalized step size allowed in the trace. Usually the user does not
need to adjust this parameter. Increasing MINDL will reduce the smoothness of the trace by
overriding the angle criteria, resulting in more aggressive projection and fewer simulation
points. Reducing MINDL will enhance the smoothness and increase the number of points in
the trace. (Default: 0.1)
• FREQUENCY specifies how often the binary output (solution) files of the trace are saved. All
I-V points are saved in outputfile. However, the PISCES solution files corresponding to
these points are saved only if they are designated by FREQUENCY. If FREQUENCY=0,
none of the solutions is saved, except perhaps the turning points (see below). If
FREQUENCY=5, e.g., the solution file of every fifth point will be saved to files named soln.5,
soln.10, etc., along with its PISCES input file (input.5, input.10,...) and output I-V file (iv.5,
iv.10,...). (Default: 0)
• TURNINGPOINTS is either YES or NO. If it is YES, the binary output (solution) file from
PISCES will be saved whenever the slope of the I-V curve changes sign, i.e., there is a turning
point. The name of the output file is soln.num,where num is the number of the current
solution. For example, if the 25th point has a different sign than the 24th point, Tracer will
save a file called soln.25. (Default: NO)
• VERBOSE is either YES or NO. If it is YES, certain information about each solution (which
the user may not be interested in) is printed in outputfile. The information consists of the
external control-electrode voltage, the load resistance on the control electrode, the slope
(differential resistance) of the solution, the normalized projected distance of the next
simulation I-V point, and the normalized angle difference between the last two simulation
points. (Default: NO)
9.3.4 Examples
– 14
1. Step-size control will begin when the control electrode’s current exceeds 1.0 ×10 A/µm. In the
input deck itlim has been set to 12. Only essential information is saved in outputfile. The
solution file of every tenth point, as well as any turning points, will be saved.
2. In a simulation with one or two open contacts, we want to keep the current through the open
– 16
electrodes below 1.0 ×10 A/µm, regardless of the current through the control electrode. Thus
RELMAX is set to a very low value so that it will not be a factor in determining the current at the
open contact(s).
9.4.2 Syntax
FIRSTSOLUTION=<char> [OPENCONT1=<int>] [OPENCONT2=<int>]
[SIMULATOR=<char>] [VOPEN1=<real>] [VOPEN2=<real>]
9.4.3 Parameters
• FIRSTSOLUTION is either INITIAL, LOAD, or CURRLOAD. In all cases a solve
statement should be present in the PISCES input deck (inputfile). The parameters of this
solve card in inputfile are not used but rather the card itself is used to mark where a
PISCES solve card should be placed by Tracer in inputfile (see SECTION 10).
Tracer run in which voltage bias conditions were used on the zero-current electrodes for the
last simulation point.
• OPENCONT1 and OPENCONT2 are the numbers of electrodes (between 1 and 9, inclusive)
with a zero-current bias condition. There can be either zero, one, or two open contacts. When a
device has an open contact, the user does not have to worry about convergence at low device-
current levels. Tracer will automatically adapt the bias conditions to guarantee convergence.
(Default: none)
• VOPEN1 and VOPEN2 must be used if and only if there is an open contact and
FIRSTSOLUTION=LOAD (voltage bias condition on open contact(s)). The values of
VOPEN1 and VOPEN2 are the voltages of the open contacts OPENCONT1 and
OPENCONT2, respectively, in the loaded solution file designated on the load card of
inputfile. If there is only one open contact, VOPEN2 should not be specified. (Defaults:
0.0)
9.4.4 Examples
1. The trace starts by solving an initial solution at zero bias and uses PISCES-2ET as the simulator.
Electrode 2 is an open contact.
2. The trace starts with a previous solution using only voltage bias conditions. In this loaded solution
the open contacts 2 and 4 have voltages of 0.641V and 0.509V, respectively.
The input deck used by Tracer, inputfile, is a standard PISCES file, but Tracer has certain
requirements. For understanding the basic flow of an input deck, consult the PISCES manual. The
mesh, region, electrode, doping, and model cards must already be present in the input deck.
Additionally, the Newton solution method must be specified in the symbolic card. Other requirements
are described below.
inputfile for each such electrode, regardless of whether Tracer is to begin with an initial solution
or a loaded solution.
Even if no contact cards are required in inputfile, a line starting with “$contact” must be
present so that Tracer will know where to add a contact statement. This contact card is necessary
because this is where the load resistance of the control electrode is specified by Tracer. There is no
problem with placing a contact card for the control electrode in the input deck as long as it does not
specify a resistance value (which should never happen). Note that at least the first five letters of
“contact” must appear for Tracer (and PISCES) to recognize it.
Another option must be specified in the method card if TMA’s MEDICI is used. In MEDICI,
if a solution is aborted MEDICI will try to solve for an intermediate solution and then retry the original
solution. This is not desirable when using Tracer since Tracer needs to keep track of aborted
solutions. Thus, “stack=0” should be specified in the method card of MEDICI so that it does not
attempt intermediate solutions. Analogously, the “trap” option should not be specified on the method
card in a PISCES-2ET deck.
Values in the next columns depend on which data are recorded. If requested in the FIXED
statements of tracefile, current values of fixed-voltage electrodes and voltage values of fixed-
current electrodes will be recorded for each solution point in outputfile. The order from left to
right is from low to high electrode number.
After the electrode information is recorded, further columns contain information about each
solution if VERBOSE=YES in the SOLVE card of tracefile. These columns are, from left to
right, external control-electrode voltage, load resistance on the control electrode, differential
resistance, normalized distance of the next projection, and the angle difference between the current and
previous solution points (see [1] for a description of these parameters).
The FREQUENCY and TURNINGPOINTS parameters in the OPTION card allow data to
be saved for certain specified solutions. In outputfile, those points which are saved are marked
with an asterisk next to the solution number. The files saved are the input deck, input.i; the I-V data
file, iv.i; and the solution file, soln.i; where i is the number of the solution in outputfile.
As of January 1994, Tracer works with PISCES-2ET, some in-house versions of Stanford PISCES,
and to some extent md3200 or md10000, MEDICI Version 1.2.2.1 Use of MEDICI is not yet robust
and thus Tracer may or may not complete a trace using MEDICI. Also, MEDICI cannot yet be used
for simulations with open contacts. Tracer has not been thoroughly tested for simulations with two
open contacts.
The capability to calculate admittances using the difference method must be added to Tracer
if it is to use simulators which cannot perform AC analysis when a load resistor is present (a previous
version has this capability, so it should not be hard to implement).
A tar file has been created which contains the Tracer source code; an executable version of
PISCES; this document in FrameMaker, MIF, and postscript formats; and a test suite of verified
examples. This test suite contains all the files needed to run the examples in the appendix.
1. These implementations were developed in connection with Advanced Micro Devices, where TMA
software is used, as part of a summer internship.
In each of the Tracer examples below, a description of the simulation is given along with the
command line used to invoke Tracer and figures with the listings of inputfile (the PISCES input
deck), tracefile, and outputfile.
13.1 BVCEO
The BVCEO experiment is conducted by biasing an npn bipolar transistor’s collector positively with
respect to the emitter while the base is left open. The PISCES input deck, bvceo.pis, shown in Figure
13.1, defines the mesh, region, electrodes, doping, emitter contact, physical models, and solution
method. Even though the contact card is not for the collector, which will be the control electrode, the
presence of the card ensures that Tracer will be able to find the correct place to insert a contact card
for the collector when it needs to. If we did not wish to use the contact card in bvceo.pis, we would
still have to insert a line beginning with “$contact” above the model and symbolic cards. Notice that
“nowarn” and “curvetrace” are specified on the options card and “newton” is specified on the symbolic
card, while nothing is specified on the solve card.
In the trace file bvceo.tra (Figure 13.2), the FIXED card sets the voltage on the emitter
electrode (num=1, as defined by bvceo.pis) to a constant value of 0.0V and states that the current
through this electrode will not be recorded in outputfile. Electrode 3, the collector electrode, is
designated as the control electrode. The CONTROL card states that the first solution will have a
solve
end
Figure 13.1 The input file bvceo.pis for the BVCEO example.
Figure 13.2 The trace file bvceo.tra for the BVCEO example.
collector voltage of 0.0V, while the second solution will have a collector voltage of 0.1V. Tracing will
continue until the collector voltage equals or exceeds 20V. If the initial step of 0.1V proves to be too
large for convergence, Tracer will cut the step size in half, possible more than once, until it converges
on a solution, and then will proceed from this solution.
In the SOLVE card, we specify that the base electrode (num=2) is to be treated as an open
contact during the trace. Also, tracing will begin with a thermal-equilibrium solution and PISCES-
2ET will be used for the simulation. Finally, the OPTION card specifies that only essential I-V data
will be saved in the output file; the PISCES iteration limit is set to 15, agreeing with the PISCES deck
in the input file; PISCES solutions will be saved for any turning points as well as for every fifth
solution point; smoothness of the I-V curve will not be enforced until the collector current is greater
– 12
than 5 ×10 A/µm; and while voltage biasing is used on the open base contact, a solution will be
– 19
accepted only if the current through the base is less than 5 ×10 A/µm (unless the RELMAX
condition predominates).
While Tracer is running, the output of the PISCES runs are sent to the standard output,
along with messages announcing when solutions are written to the output file. The output file, named
bvceo.out in the command line, is shown in Figure 13.3, and a plot of the collector current vs. collector
voltage is shown in Figure 13.4. In bvceo.out, we see that every fifth solution, along with solutions 24
and 36 (the turning points), has been saved in files named soln.5, soln.10, etc. Additionally, the last
solution was saved in the file soln.last, although there is no asterisk marking the last solution in
bvceo.out.
Figure 13.3 The output file bvceo.out for the BVCEO example.
Figure 13.4 Collector current vs. collector voltage for the BVCEO example.
At the top of bvceo.out, column headings mark the solution number, control-electrode
(collector) voltage, control-electrode current, open-contact (base) voltage, and open-contact current as
Soln, Vctrl, Ictrl, Vcurr, and Icurr, respectively. We see that the collector voltages for the first, second,
and last solutions are 0.0, 0.1, and 20.18V, respectively. The final solution does not have a collector
voltage of exactly 20V, as specified in bvceo.tra, because Tracer only guarantees that the curve will
be traced out to at least 20V, not exactly 20V.
Other information regarding the trace must be inferred from the PISCES output displayed
while Tracer is running (not shown). From this output we can see that voltage biasing was used on
the open base contact for the first few solutions, in which the collector current is too small to allow
stable use of zero-current biasing. A few PISCES simulations are actually run for each I-V point, with
minor adjustments on the base voltage being made until the base current is less than ABSMAX. When
the collector current is large enough, Tracer places a zero-current bias on the base. We can also see
that a variable load resistor is placed on the collector when the collector current exceeds MINCUR.
After this, the step sizes are regulated to produce a smooth curve.
For Tracer, another PISCES input deck must be created to use as the input file (Figure
13.6). In mesvg.5.pis the mesh file generated by mes.pis, mes.mesh, is read in, preempting the mesh,
eliminate, region, electrode, and doping cards. Since Tracer will be starting with a previous solution,
the name of the solution file to load must be given in mesvg.5.pis. This load statement appears directly
above the solve card with the file name mesvg.5.ini, the solution file generated by mes.pis.
The trace file mesvg.5.tra is shown in Figure 13.7. In the three FIXED cards, the voltages of
the source and substrate (num=1 and num=4, respectively, as defined by mes.pis) have been fixed at
0V, while the gate voltage (num=2) has been fixed at -0.5V. The current through the gate electrode
will be recorded for each solution in the output file. The CONTROL card of mesvg.5.tra specifies that
the drain (num=3) will be swept from 0.0V to a voltage where the current is greater than or equal to
–4
4.1 ×10 A/µm, with an initial drain voltage step of 0.2V. On the SOLVE card, FIRSTSOLUTION
is specified as LOAD, consistent with the input file mesvg.5.pis, and PISCES-2ET is designated as
the simulator to use. Since VERBOSE is NO on the OPTION card, only the essential I-V data will be
recorded in the output file. The iteration limit is 30, consistent with mesvg.5.pis, and every ninth
solution, as well as those corresponding to turning points, will have its solution file saved.
title mes.pis
$*** regions
region num=1 ix.lo=1 ix.hi=53 iy.lo=4 iy.hi=41 gaas
region num=2 ix.lo=1 ix.hi=53 iy.lo=1 iy.hi=4 oxide
region num=2 ix.lo=16 ix.hi=34 iy.lo=1 iy.hi=7 oxide
$*** doping
dop ascii x.l=0.0 x.r=10.0 inf=mei.dop
dop gaus x.l=-1 x.r=3.0 dos=5.0e13 cha=0.0607 peak=-0.0709
+ n.t erfc.lat lat.cha=0.0866
dop gaus x.l=7.0 x.r=11 dos=5.0e13 cha=0.0607 peak=-0.0709
+ n.t erfc.lat lat.cha=0.0866
$*** material
material num=1 eg300=1.42 affinity=4.07 vsat=10.0e6
+ permi=13.1 nc300=4.35e17 nv300=8.35e18
interface qf=-1e12 x.min=0.0 x.max=10 y.min=-.01 y.max=6.0
$*** contact
contact num=2 alu workf=4.84 surf
solve ini
symb newton carrier=2
solve v2=-0.25
solve v2=-0.5 proj outfil=mesvg.5.ini
end
title mesvg.5.pis
mesh inf=mes.mesh
load infil=mesvg.5.ini
solve
end
Figure 13.6 The input file mesvg.5.pis for the GaAs MESFET example.
Figure 13.7 The trace file mesvg.5.tra for the GaAs MESFET example.
Figure 13.9 shows the output file, mesvg.5.out, in which the solution number, drain voltage,
drain current, and gate current have been recorded as Soln, Vctrl, Ictrl, and I2, respectively. The
solution files of points 9, 18, 27, and 29 (a turning point), as well as of the last point (not marked in the
output file) were saved as soln.9, soln.18, soln.27, soln.29, and soln.last, respectively. A plot of the
drain current vs. drain voltage is shown in Figure 13.8.
Figure 13.8 Drain current vs. drain voltage for the GaAs MESFET
Figure 13.9 The output file mesvg.5.out for the GaAs MESFET.
Reference
[1] R.J.G. Goossens, S.G. Beebe, Z. Yu, and R.W. Dutton, “An Automatic Biasing Scheme for
Tracing Arbitrarily Shaped I-V Curves,” IEEE Trans. on Computer-Aided Design, vol. 41, pp.
310-317, March 1994.
The prototype to the mixed-mode simulator was first developed in the summer of 1992 by Hui Wang,
a visiting scholar from Tsinghua University, Beijing, China [1]. Funding for this project came from
AASERT contract DAAH04-93-G-0178 on the parent Army contract DAAL03-91-G-0152. The
MESFET example is provided by Matsushita Electric Industrial Co. in Japan and the GaAs/AlGaAs
LED example is from OCD, HP, in San Jose, CA.
14.1 Introduction
The Stanford mixed-mode interface provides a mechanism to include complex numerical devices in
the SPICE circuit simulator when compact models are either inadequate or not available. Such devices
include GaAs MESFET’s, heterostructure devices, short channel MOSFET’s, and optical devices.
SPICE3f2 developed at UC Berkeley is used as the circuit simulator upon which the numerical device
model is added.
The mixed-mode simulator comes configured to work with Stanford’s version of PISCES-2ET.
However, the modularity of the interface between SPICE and PISCES is such that any other device
simulator may be interfaced with SPICE. Refer to “System Reconfiguration for a Different Device
Simulator” on page 303 for more information.
SPICE has many different types of analyses which Stanford’s mixed-mode simulator utilizes.
Currently, a numerical device may be used in DC, AC small-signal, and transient analyses. Pole zero
analysis has been included; however, this feature tends to have trouble in convergence for all but the
simplest of circuits. Sensitivity analysis will be added in a later release.
F(V ) = 0 (14.1)
where V is the vector of node voltages and F represents the sum of the currents into each node in the
circuit, and both V and F have dimension of N, the number of nodes in the circuit.
Applying Newton-Raphson to the above equation yields the linear matrices shown in Equation 14.2
where J ( V ) is the Jacobian as given in Equation 14.3.
V i + 1 = V i – J –1 ( V i )F ( V i ) (14.2)
∂F 1 ∂F 1 ∂F 1
…
∂V 1 ∂V 2 ∂V N
J (V ) = … … (14.3)
∂F N ∂F N ∂F N
…
∂V 1 ∂V 2 ∂V N
i i+1
For each Newton iteration, the previous voltage V is known and hence V can be computed. A
circuit interpretation of Equation 14.2 and the Jacobian matrix (Equation 14.3) can be made as follows.
Because each Jacobian element has units of the conductance, we hereafter use G to replace J. By
multiplying G on both sides of Equation 14.2 from the left, one obtains the following set of linear
equations at ( i + 1 ) -th iteration where i starts from 0:
i i i i
G Vi + 1 = G V – F (14.4)
i i i
The matrix G and vector F are evaluated at V . The above equation represents a linear circuit as far
i+1 i
as node voltages V because both the coefficient G and the source term, i.e., the right hand side
i+1 i
(RHS) of the above equation, are independent from V . Furthermore, G can be interpreted as the
linear conductance components, including both the linear components in the original circuit and
i i i i
differential conductance at V , in the equivalent circuit. G V – F can be viewed as the current
sources. For a more detailed discussion, refer to McCalla [2].
The NR iteration is terminated when the convergence is reached, i.e., the change in V between two
consecutive iterations is smaller than a predefined tolerance. In SPICE, it is required that the current
change in each circuit branch is also below certain criterion when the convergence is considered to be
reached.
Transient analysis is performed in a similar manner. For each time step, Newton-Raphson iterations
are performed until convergence is met. In addition, the truncation error due to the time discretization
is checked to determine if the time step is acceptable in terms of accuracy. If this error is too large, the
time step is reduced and the computation is repeated.
For AC analysis, the DC solution is first sought and then the non-linear devices are linearized at the
solution. This small signal equivalent circuit is used to find the AC response given an AC excitation
vector.
The device simulator is responsible for supplying the terminal currents and the admittance matrix at
the given voltage boundary conditions. The device simulator first solves F ( w, V , t ) = 0 where
w = f ( ψ, n, p ) for the drift/diffusion model and w = f ( ψ, n, p, T n, T p, T L ) for the duel energy
transport model (see PART I) V is the terminal voltage on the device. From the solution, the terminal
current density and conductance matrix can be calculated. For AC analysis, the admittance matrix must
be calculated using frequency dependant AC analysis which is available in PISCES.
This interface utilizes a two level Newton scheme to solve a circuit that includes numerical devices.
On the upper level is the circuit iterations to determine the node voltages. For each one of these
Interface
PISCES
Device Simulator
I = solution to semiconductor
equations and models
PISCES-IIB
SPICE Write
PISCES-2ET
Numerical Device
Simulator Vendor PISCES
Find Bias for Input Deck
Numerical
Devices Hal’s Simulator
Sockets
SPICE-TCAD
Execute Interface Execute TCP
Numerical Device Numerical Device
Runs Simulation PVM
TCL
Load Device
Simulation Results in PISCES-IIB
G Matrix and I Vector
Read PISCES-2ET
Numerical Device
Simulator Results Vendor PISCES
Hal’s Simulator
iterations, the device simulator solves for the terminal currents and conductance matrix using another
Newton iteration within the device simulator (the lower level). The advantage of this method is the
modularity created in the code.
Figure 14.2 shows the flow diagram for the execution of a single iteration step at the upper level. The
blocks on the far right hand side represents different modules for the device simulator and/or method
of execution. The sequence of steps is as follows.
1. Determine the voltage boundary conditions, time step, and/or frequency SPICE
requires for the current iteration/solution. The exact information depends on
the type of analysis.
2. Create the numerical device input deck that solves for the information
requested by SPICE. A separate module exists for each device simulator that is
interfaced to SPICE.
5. Based upon the analysis type, load the SPICE conductance (admittance) matrix
and current vector.
Since a device simulation is required for each circuit iteration, the interface algorithms are given some
intelligence to limit the computations required by the numerical device. The most important aspects or
those algorithms are as follows:
In the full level Newton algorithm, the circuit matrix and device matrix are combined into a single
matrix and all variables are solved simultaneously. As a result, the computational time is reduced, but
there are trade-offs.
1. It has better convergence for DC analysis if node voltages are unknown. The
full Newton requires all circuit node voltages to be specified within a certain
percentage; otherwise, it fails to converge.
2. The modularity in two level Newton allows many device simulators to be used
simultaneously. For example, PISCES may be used for standard devices in the
circuit whereas an in house device simulator may be used for novel devices
Likewise, there are a number of disadvantages to the two level Newton algorithm:
1. Given a good estimate of all the circuit node voltages, the full Newton
algorithm converges more quickly than the two level Newton.
2. Similarly, since transient analysis involves small voltage changes at each time
step making the problem well behaved, the full Newton method converges
much more quickly for this case as well. Mayaram determined a factor of 1.7
times as fast [3].
14.6 Parallelization
In order to reduce the overall user time per circuit iteration, the numerical devices (i.e. the lower level
Newton iterations) can be solve on different nodes on a network of machines. Because each device
simulation is an independent unit, parallelization is relatively trivial.
The mixed-mode interface contains two mechanisms by which the parallel simulations are executed.
The most basic of these algorithms utilizes standard Unix sockets. A more elegant mechanism involves
PVM (Parallel Virtual Machine) which allows a heterogenous network of computers to act like a
parallel machine [5]. Figure 14.2 shows a block diagram of how PVM is utilized.
The mixed-mode interface spawns a queueing process on a node of the virtual machine. This process
ascertains the number of nodes that can be used for a device simulation. Whenever the interface
program needs solutions from a number of numerical devices, it sends a request to the queueing
Node 1
Node 2
Queueing
Mixed-Mode Controller
SPICE Interface for PVM Node 3
Execution
Node 4
Sequential
Execution
Node 5 Parallel
Virtual
Machine
program which then executes the device simulations on the available nodes. Once the device
simulations are complete, the queueing program sends back the solutions to the interface. As a result,
the computational load on the nodes is optimized.
15.1 Introduction
The Stanford mixed-mode simulator provides a mechanism to incorporate numerical devices into a full
circuit simulation with linear and non-linear components. The mixed-mode simulator uses the industry
standard SPICE as the circuit simulator with an added model for the numerical device. This numerical
model is configured modularly such that any device simulator may be added. (Refer to “System
Reconfiguration for a Different Device Simulator” on page 303 for more information.) For the purpose
of this document, Stanford PISCES is used.
This chapter describes how to include a numerical model in a circuit simulation. All descriptions are
based upon SPICE version 3f and hereafter referred to as SPICE only.
As with all element cards in SPICE, this card describes a specific instance
of the numerical device in the circuit. The name of this instance is used as
the predecessor for all files accessed in solving the specific numerical
device.
SYNTAX
PARAMETERS
Nxxxxxxxx
The N identifies the element type as being a numerical device and the x’s
are used to represent a unique character string, which can be up to eight
characters long, to identify the device.
n1 n2 . . . nM
The node numbers in the circuit to which the nM nodes of the numerical
device are connected. They correspond in sequence to the electrodes of the
device as specified in PISCES. There is a maximum limit of ten nodes per
device. The voltage on any node is specified with respect to the last node as
given by the nodes parameter on the .model card.
mname
The model name links this instance of the numerical device to the
appropriate .model card.
area
The area or length factor for the device. All currents, capacitances, and
conductances are multiplied by this value. The default is 1.0.
off
If off is specified, the DC operating point is determined with the terminal
voltages for that instance of the device set to zero.
ic
These initial conditions for the device are used with the UIC specification
on the .tran card. Each voltage initial condition is specified with respect to
the last node.
temp
Operating temperature of the device. If none is given, the circuit
temperature is used as the default.
EXAMPLES
Nbjt1 8 7 5 npndev 5
This adds a numerical device with the name of Nbjt1 connected to circuit
nodes 8, 7, and 5. These nodes correspond to electrode numbers 1, 2, and 3
in the numerical device’s mesh. Electrode number 3 is taken as ground.
Analogously with SPICE’s compact model for a bipolar junction transistor,
Electrode 1 (Node 8) is the collector, Electrode 2 (Node 7) is the base, and
Electrode 3 (Node 5) is the emitter. An area factor of 5.0 scales the device
appropriately.
BUGS
Some software paths have not been tested and may result in problems.
As with all .model cards in SPICE, this one is used to describe all
instances of the same type. Common parameters are specified and used by
all instances that reference the given model name.
SYNTAX
PARAMETERS
mname
Model name used for reference on the element cards. The mixed-mode
interface expects to find certain files with this name as the predecessor to
the file names.
mtype
The model type is given by a word mnemonic from Table 15.1. The
mnemonic indicates the polarity of the numerical device and the
mechanism by which the numerical simulation is executed.
For the pvm model, the mixed-mode interface uses PVM message passing
to send the biases to a queueing routine and ultimately a “slave” routine on
a node of the parallel virtual machine. Refer to “Running a Mixed-Mode
Simulation” on page 286 for more information.
Table 15.1
mtype Description
pis generic numerical device (default)
npis n polarity numerical device (npn, nmos, pn junction)
ppis p polarity numerical device (pnp, pmos, np junction)
prl generic numerical device simulated in parallel
nprl n polarity numerical device simulated in parallel
pprl p polarity numerical device simulated in parallel
pvm generic numerical device simulated using pvm
npvm n polarity numerical device simulated using pvm
ppvm p polarity numerical device simulated using pvm
nodes
The number of numerical device nodes. This parameter defaults to 10 if no
value is given. SPICE does not have a problem if this number is larger than
the actual number of electrodes on the numerical device. The unidentified
nodes are connected to ground by default. Likewise, one must be certain
that the numerical device simulator can handle a “solve” statement with
voltages specified for electrodes that do not exist. To avoid this problem,
one should always specify this parameter.
The nodes parameter may be less than the actual number of nodes in the
numerical device. In this case, the extra nodes on the numerical device are
considered floating and always ignored by the SPICE circuit simulator.
Other .model parameters related to floating nodes include float, vbc, and
bcdep. These parameters allow the user to switch boundary conditions on
the floating node during a DC circuit solution.
dtype
This integer value provides SPICE with some a-priori knowledge of the
device behavior. This parameter is not necessary, but is useful in achieving
circuit convergence. Specifying a device type allows SPICE to make good
guesses and limit voltage changes across pn junctions. The supported
devices and their integer representations are configured in the NPISCdefs.h
file. The default configuration is given in Table 15.2. The mtype parameter
determines the polarity of the device. The electrode sequence for the
electrodes is analogous to the SPICE compact model equivalents.
Table 15.2
vt0 or vto
This parameter represents the absolute value of the turn on voltage for the
device. For a bipolar junction transistor it is Vbeon, for a diode it is Vdon,
and for a FET it is Vgson. The default value at room temperature is 0.7 V
which is typical for a silicon pn junction. The polarity of the device model
determines the sign.
vmax
This is the maximum voltage (in signed real number) placed on any
numerical device electrode with respect to the reference electrode. The
default value is 50.0 volts which one needs to increase for power devices.
vmin
This is the minium voltage (in signed real number) placed on any numerical
device electrode with respect to the reference electrode of that numerical
device. The default value is -50.0 volts.
hifreq
This value separates low and high frequency simulation of the numerical
device. When the circuit frequency is below hifreq, the device’s small
signal parameters are calculated from a zero frequency device simulation.
When the circuit frequency is above hifreq, the device’s small signal
parameters are calculated from a small signal device simulation at the given
frequency. The trade-off is accuracy versus simulation time. The default
value for hifreq is 1 kHz.
maxtrys
The maximum number of attempts the mixed-mode interface tries to get the
device simulation to converge. Different bias stepping schemes are
attempted with each try.
level
This integer represents a specific device simulator as defined in the
NPISCdefs.h. Currently, the only one supported (and the default) is
Stanford PISCES which is defined as level=0. This integer value
determines how the method and model parameters are interpreted as well
as the calling of the correct device simulator. More levels may be added
model
This is an integer value for the binary encoded representation of the model
card in the device simulator input deck. The meaning of the integer value is
described in the NUMDEVdefs.h file and interpreted by a routine
specifically written for that device simulator. If this parameter is not
specified and the file mname.model exists, the model card is taken as
the content of that file. If neither is given, this card is not included in the
input deck for the device simulator. Refer to “Running a Mixed-Mode
Simulation” on page 286 for more information.
method
This is an integer value for the binary encoded representation of the method
card in the device simulator deck. The meaning of the integer value is
described in the NUMDEVdefs.h file and interpreted by a routine
specifically written for that device simulator. If this parameter is not
specified and the file mname.method exists, the method card is taken
as the content of that file. If neither is given, this card is not included in the
input deck for the device simulator. Refer to “Running a Mixed-Mode
Simulation” on page 286 for more information.
currents, one uses a current boundary condition on the floating node. The
value of the voltage and current source on the floating node is zero.
The parameter float contains the integer number for the floating node. It
must be greater than the number specified by the nodes parameter. The
parameter vbc determine the voltage at which the boundary condition is
switched. The parameter bcdep specifies the positive and negative
electrode on which the boundary condition depends. Refer to the examples
for clarification.
oedev
An integer value of one for oedev tells SPICE that the numerical device
will produce optical output. As a result, SPICE adds an equation to its
circuit matrix to store the value of that output. In a future release, an integer
value of two will tell SPICE that the numerical device has light incident
upon it. The amount of incident light will be determined by referencing a
circuit node. The default of zero means no photons are involved.
tnom
This parameter specifies the nominal temperature at which the parameters
are given. Currently, the only parameter affected by tnom is vt0. The
default value is 300K.
EXAMPLES
parameter. The turn-on voltage for the device is 0.7 volts. The model and
method for the device simulation is specified on the card and is explained
in the next section. The number of nodes is given as three which
corresponds to the three nodes on the previous element card.
In addition, this device will produce optical output and hence, an optical
response. The model and method lines for the numerical simulator input
deck is given in the files hpled.method and hpled.model.
BUGS
The advanced parameters for floating layers have been generalized, but
may not be useful for all cases.
When vmax and vmin are set stringently, they tend to cause SPICE to
oscillate between the extremes. These parameters should be set such that
they prevent the numerical device from entering breakdown.
Some software paths have not been tested and may result in problems.
mynetlist.spi
This file contains the circuit netlist of all the element cards and model cards for the circuit. Any valid
SPICE construct is permitted. The name of the file can be arbitrary, but it is recommended that the
extension .spi is used so that it may easily be recognized.
mynetlist.raw
The easiest way to run SPICE is to create a net list such that the circuit simulation is performed in batch
mode. As a result, the user needs to save the solution SPICE computes. The name of the file can be
arbitrary, but it is recommended that the extension .raw is used. To run SPICE in batch mode one
types the following on the command line:
mname.msh or mname.mesh.pis
One of these files must be supplied to load or create the mesh during the PISCES simulation. The first
is a standard PISCES mesh file. The second contains valid PISCES commands for creating the mesh
from scratch. The mname corresponds to the model name given on the SPICE .model card. The
extension .msh or .mesh.pis must be exactly as shown. The mixed-mode code first looks for
mname.msh, if that file does not exist, it looks for mname.mesh.pis. If neither exists, SPICE
aborts.
Although these files are based on PISCES, they should be completely compatible and/or analogous for
other device simulators. The actual format is determined by the programmer who adds the new device
simulator.
Using mname.msh provides for faster simulation because PISCES just loads the mesh and doesn’t
have to create it each time. However, for some simulations, one may choose to create the mesh each
time and therefore, the second option is available. When using mname.mesh.pis, the user should
not include an end card.
mname.soln.init
This is a PISCES solution file containing the initial zero bias solution for the model mname. This only
has to be found once for the model since all instances use the same starting point.
Nxxxxxxxx.pis
This file contains the PISCES input deck that the mixed-mode interface writes for the instance named
Nxxxxxxxx.
Nxxxxxxxx.out
This file contains the output of the last PISCES simulation run for the specific instance.
solution for iteration before the previous. For transient analysis, .prev0 file contains the solution for
the current time step, the .prev1 file contains the solution for the previous time step, and the .prev2
file contains the solution for the time step before the previous time step.
+1 turns on srh
+2 turns on consrh
+4 turns on auger
+8 turns on bgn
+16 turns on conmob
+32 turns on analytic
+64 turns on fldmob
+128 turns on surfmob
+256 turns on impact
+512 turns on ccsmob
+1024 turns on fermi
+2048 turns on incomplete
+4096 turns on photogen
For example, to generate the following model and method lines in the PISCES input deck,
1. 2nd and tauto must be turned off for mixed-mode simulations using PISCES.
The hostname.prl file contains a listing of the hosts on which the device simulations are executed.
Each numerical device in the circuit is assigned a machine from the list and it is always executed on
that machine. If there are more devices than machines then some machines can have more than one
device executed on it. It is recommended that the machines be listed in decreasing computational
power.
The login and password for the remote machine is required. The mixed-mode code looks for the host
name in the .netrc file. However, since the login and password for many hosts are the same, the user
can use a machine name of “spice” in the .netrc file and the mixed-mode code uses that login and
password for all hosts not found in the .netrc.
Three files are necessary for using PVM to execute a mixed-mode simulation. The user should place
the executable files npiscctrl and npiscslave in the PVM bin directories of the appropriate
machine architectures. In addition, hostname.pvm needs to contain a listing of all the nodes used
for the mixed-mode simulation; however, they do not need to be added to the parallel virtual machine.
If they are not in the parallel virtual machine, the mixed-mode simulator attempts to add them.
Finally, the user should start and/or make sure that the PVM daemon is running on the local machine
and then start the SPICE simulation. If the PVM daemon is not running, mixed-mode defaults to
sequential execution.
16.1 Introduction
This chapter describes examples using three types of analyses. The first is a CMOS inverter analyzed
with a DC sweep on the input. The second is a GaAs MESFET analyzed for its high frequency
response. The third is a six transistor SRAM cell analyzed through a read/write cycle. Another
transient analysis is used to analyze a fiber optic transmitting circuit which includes a heterostructure
LED with a floating layer.
For this simulation the mesh files for the two devices are created separately and are called
fmrpmos.msh and fmrnmos.msh. These mesh files contain the grid and doping for the
numerical device. The PISCES method and model options are chosen by the binary encoded SPICE
parameters method and model as follows.
vin 4 0 dc 0
vdd 1 0 dc 5
For the DC analysis, the input voltage is swept from zero volts to five volts and the output voltage and
supply current are monitored. The output waveform is as what would be expected from a CMOS
inverter and is not shown.
For this simulation, the file mes.mesh.pis contains the PISCES cards for creating the MESFET’s
mesh during each device simulation. Note that nothing specific has to be written in the SPICE netlist
Igate Idrain
MESFET
+ +
Vgate Vdrain
- -
* MESFET circuit
vgg 9 0 dc 0.0 ac 1
vdd 10 0 dc 2.0
vgate 9 8 0
vdrain 10 1 0
ld 1 2 0.08n
rd 2 3 1.52
cd 2 0 212f
lg 8 7 0.105n
rg 7 4 0.92
cg 7 0 208f
ls 0 6 0.005n
rs 6 5 1.38
Figure 16.2 Two port network and netlist for AC analysis to find the
Y-parameters
to differentiate using a previously created mesh file. In addition, the files mes.model and
mes.method are also supplied and each contains one of the following lines respectively.
* S11
x S22
Figure 16.3 S-parameters for GaAs MESFET in the frequency range of 1GHz
to 40GHz.
An AC analysis is performed at each input of the two port network with other shorted. The
Y-parameters are calculated by taking the currents Igate and Idrain and dividing by the AC excitation
magnitude. A simple transformation yields the S-parameters for a 50 ohm cable. The results are
displayed in Figure 16.3.
The transient analysis of the circuit consists of a write and read cycle. The output of the flip-flop and
the data lines are observed. Figure 16.5 contains the net list for the circuit. There is an extra compact
NMOS transistor for leveling the charges of the data line capacitances. After a read or write cycle, this
VDD
L
C L
C
R
WB
WB
CS
transistor turns on to allow the charge between the data lines to distribute evenly. The data lines must
be equal before anymore actions can take place.
The models for the numerical devices are the same as those used for the CMOS inverter. The compact
model is SPICE Level 2 model which closely represents the numerical model. The devices in the
circuit are not the state-of-the-art, but they do provide a good example. The SPICE .nodeset line is
* supply line
Vdd 1 0 5
* pass gates
NM5 4 7 2 0 fmrnmos area=2
NM6 5 7 3 0 fmrnmos area=2
* control lines
M7 1 1 2 0 ndev l=2u w=25u
M8 1 1 3 0 ndev l=2u w=25u
M9 2 9 6 0 ndev l=2u w=25u
M10 3 8 6 0 ndev l=2u w=25u
* access line
M11 6 10 0 0 ndev l=2u w=50u
* level line
M12 2 12 3 0 ndev l=2u w=25u
* compact models
.model pdev pmos
+ level = 2.000 vto = -0.5677 gamma = 0.51
+ phi = 0.73 tox = 2.5000E-08 nsub = 4.0251E+14
+ nfs = 6.8412E+11 xj = 2.5976E-8 ld = 1.8553E-07
+ uo = 1.705E+02 ucrit = 2.2295E+05 uexp = 0.2918
+ vmax = 6.8950E+04 neff = 9.285E+1 delta = 0.5587
+ cj = 2.7923E-04 cjsw = 2.1805E-10 pb = 0.7316
+ mj = 0.4729 mjsw = 0.2195 fc = 0.5
+ cgdo = 4.71e-10 cgso = 4.71e-10
* numerical devices
.model fmrpmos npvm dtype=2 vt0=0.8 model=94 method=336 nodes=4
.model fmrnmos npvm dtype=2 vt0=0.8 model=94 method=336 nodes=4
used to set the flip-flop to an initial state rather than having it start indeterminately. The .tran card
specifies the limits of the transient analysis.
The numerical device models are specified as pvm. As a result, each numerical device simulation is
executed on a different node of a parallel virtual computer. The file hostname.pvm is provided
with a list of hostnames for the nodes. The network uses a common file server whose disks are mounted
on each of the remote hosts.
5
4 C
Data Lines
3
2
C
1
0
L
5
Latch
4
3
2
1 L
0
Read/Write Cycle
5
4 Write 1 Read 1
3
2
1
0
0 5 10 15 20 25 30 35 40
Time (ns)
Figure 16.6 shows the output waveforms of the SRAM cell. It is important to note how fast the
difference between the data lines becomes small. This time determines the minimum read/ write cycle
time.
cathode
contact emitted light
n
n
p
floating
layer n+ p
anode
contact
Figure 16.7 Cross section of LED structure that is used in a fiber optic
transmitter circuit.
At low currents and biases, the voltage boundary condition introduces some error, but it is tolerable
since that is not the range of interest for the operation of the device. For higher currents, the node is
kept floating by placing a current boundary condition with a current source value of zero. In order to
handle this switching, there are three special parameters.
This device is placed into the transmitter circuit shown in Figure 16.8. The purpose of this circuit is to
take an ECL signal, convert to TTL, and either turn on or off the LED. The SPICE netlist for the circuit
is also given.
First, one notes that this device is called using ppis model because a ppis device is defined such that
the cathode is electrode 1 and the anode is electrode 2. An npis device has the reverse. In general, the
n-polarity devices are such that the device is in its SPICE standard forward operating state.
Vcc
+5 V
10 µF
0.1 µF
33 Ω 33 Ω
TTL
in 56 µF 270 Ω
* supply
Vcc 6 0 DC 5
* driver circuit
R8 2 3 33
C4 2 3 75p
R9 3 4 33
r10 4 0 270
Vhp1 6 5 DC 0
nhp1 4 5 hpled
.options abstol=1e-9
.model hpled ppis vt0=1.4 maxtrys=2 nodes=2 dtype=3 float=3
+ vbc=1.0 bcdep=21
Figure 16.8 Circuit schematic and netlist for LED transmitter simulation.
The rest of the parameters on the .model line are set to account for the heterostructure device. vt0 is
set to turn on the device. The maxtrys parameter limits the maximum number of tries to two. This
parameter is set low because the cutting of the bias step does not take into account a changing boundary
condition. The number of nodes is set to the actual number of contacts in the SPICE circuit. The
floating node is not connected to anything and hence, SPICE does not need to know about its existence.
The three parameters float, vbc, and bcdep specify the node and the switching point for the boundary
condition.The floating node is electrode number 3 as given by the float parameter. The voltage
boundary condition is used when V21 is below 1.0 volts and the current boundary is used when V21 is
greater than or equal to 1.0 volts. The parameter bcdep specifies V21 as the controlling voltage and
the parameter vbc specifies the voltage of V21 when the boundary condition change is to take place.
The device uses cylindrical coordinates and hence, the area factor is allowed to default to 1.0. The
mesh is created each time the device is solved because the version of PISCES being used can not save
the heterostructure parameters. Hence, the file hpled.mesh.pis must be supplied. In addition, the
model and method cards are given in the files hpled.model and hpled.method and are as
follows.
A communications circuit engineer is most interested in the turn-on turn-off characteristics of the
diode. Therefore, the transient simulation consists of a pulse input to switch the diode on and back off.
The TTL nand gates used for supplying the current are substituted with a voltage source and a resistor.
The voltage source takes into account the finite rise and fall time of the nand gates while the resistor
simulates the output resistance of the nand gate.
Figure 16.9 shows the response of the diode. This current is directly related to the amount of photons
out for the device.
120
100
80
60
Ihp1 (mA)
40
20
0
-20
-40
-60
0 5 10 15 20 25 30 35 40
Time (ns)
Figure 16.9 Time transient response of the current flowing through LED in a
transmitter circuit.
When referring to a directory, the home directory is to be taken as spice3f/src unless otherwise
specified.
Table 17.1
inp2n.c makedefs npisc.c
npiscacload.c npiscask.c npiscconvtest.c
npiscctrl.c npiscdefs.h npiscdelete.c
npiscexec.c npiscexecprl.c npiscexecpvm.c
npiscexecseq.c npiscgethost.c npiscgetic.c
npiscitf.h npiscload.c npiscmdelete.c
npiscmethod.c npiscmodel.c npiscmparam.c
npiscmvfiles.c npiscparam.c npiscpvmmesg.c
npiscpzload.c npiscreadpvm.c npiscreadstanford.c
npiscsetup.c npiscslave.c npisctemp.c
npisctrunc.c npiscwritestanford.c numdevdefs.h
makefile.3d2
These are subject to change for later versions, but differences should not be significant.
3. In include/inpdefs.h, there are two listings for the routines that are used
to read an individual SPICE device type. The following lines are added to the
appropriate list.
void INP2N( GENERIC*, INPtables*, card* );
void INP2N();
4. The following line is added to the subroutine inp_numnodes() in lib/fte/
subckt.c. The number that is returned is the maximum number of nodes that
the numerical device can support.
case: ‘n’: return(10);
5. Add the file inp2n.c to the directory lib/inp. Modify the makedefs file
in that directory to include inp2n.c in the CFILES definition and to include
inp2n.o in the COBJS definition. Finally, add the following entry to the list
of files to make.
inp2n.o: inp2n.c
6. The following lines are added to lib/inp/inpdomod.c. These are the valid
model types available for a numerical device. The administrator can eliminate
an execution method by not including it in this listing.
} else if ( (strcmp(typename,”pis”) == 0) ||
(strcmp(typename,”npis”) == 0) ||
(strcmp(typename,”ppis”) == 0) ||
(strcmp(typename,”pvm”) == 0) ||
(strcmp(typename,”npvm”) == 0) ||
(strcmp(typename,”ppvm”) == 0) ||
(strcmp(typename,”prl”) == 0) ||
(strcmp(typename,”nprl”) == 0) ||
(strcmp(typename,”pprl”) == 0) ){
type = INPtypelook(“Npisces”);
if (type < 0) {
err=INPmkTemp(“Device type Npisces not available in this binary\n”);
}
INPmakeMod(modname,type,image);
9. In the main configuration file, one needs to include “npisc” in the listing of
devices to compile. This configuration file is located in the conf directory
which is at the same level as the src directory.
10. The user should clean everything and recompile SPICE. There may be
dependencies that are not taken into account when certain header files are
modified.
1. Modify the file lib/ckt/cktbindn.c to include the cases for nodes six
through ten. Follow the pattern already prescribed for nodes one through five.
3. The user should clean everything and recompile SPICE from scratch. There
may be dependencies that are not taken into account when certain header files
are modified.
3. Copy (or link) the files npiscslave.c and npiscctrl.c to the directory
src/bin.
4. The makedefs file in the directory src/bin has the following lines added
so that the npiscslave and npiscctrl routines are compiled and installed
correctly. In the list of CFILES, add the files npiscslave.c and
npiscctrl.c. In the list of INSTALL_TARGETS add the following lines.
$(PVM_BIN)/npiscslave \
$(PVM_BIN)/npiscctrl
A variable COBJS is defined below CFILES as follows.
COBJS = npiscslave.o npiscctrl.o
Finally, add the following line to the list of object files.
npiscslave.o: npiscslave.c
npiscctrl.o: npiscctrl.c
The first and most obvious difference is that SPICE3d capitalizes the identifier of each module. As a
result, anywhere npisc is specified, it becomes NPISC. For example npiscload.c becomes
NPISCload.c and DEV_npisc becomes DEV_NPISC.
Following are the differences between SPICE3f and SPICE3d given in the same sequence as the
original installation instructions.
1. Instead of changing the bin/config.c file, one needs to make the given
change in the bin/spice.c file.
2. Instead of changing the bin/config.c file one needs to make the change in
the include/devices.h file.
5. The file INP2N.c should be added to the lib/INP directory and the
appropriate change should be made to the Makefile in that directory.
8. The directory for the numerical device is called lib/DEV/NPISC, the files
should have “npisc” changed to “NPISC”, and makefile.3d2 is renamed to
Makefile. In addition, there are a few changes to the code as outlined here.
First, the include filenames need to be changed since SPICE3d uses a different
9. For SPICE3d, the NPISC device must be defined in the variables SUBDIRS
and MSCSUBDIRS in lib/DEV/Makefile.
10. The user should clean everything and recompile SPICE to make sure that all
dependencies are taken into account.
11. PVM has not been tested with SPICE3d and the procedure for adding the
compilation of npiscslave.c and npicsctrl.c is unknown. One can
attempt to define the appropriate variables as discussed in Section 17.1.4 on
page 306 and compile the routines manually.
Most simulators posses these attributes except for finding the low frequency conductance matrix in DC
and transient analysis. The changes to obtain this information is minimal and should not be a major
problem.
npiscdefs.h
This file is used to define the new device simulator and to give it an appropriate identification number.
This number corresponds to the level given on the .model card. All code specific to that device
simulator should be surrounded by an #ifdef structure so that it may easily be excluded in a
compilation.
npiscexec.c
The purpose of this routine is to execute the device simulations. SPICE passes this routine the voltage
boundary conditions on each electrode of the numerical device. This routine returns the solution for
the current vector and admittance matrix.
npiscexec calls two subroutines that are specific for the device simulator as determined by the level
value. One routine writes the input deck for the device simulator and the other reads the solution from
the files created by the device simulator. For Stanford PISCES, these routines are as follows and are
located in a “case” structure in the file.
npiscwriteMYSIMULATOR.c
This routine is responsible for writing the device simulator input deck. A pointer to a structure
containing the information about the bias is passed to the routine. An error code or 0 is returned from
this subroutine call.
The bias structure contains information about the SPICE instance name, SPICE model name, voltage
bias, time step, frequency, temperature, etc. (Refer to numdevdefs.h for all the information
contained within the structure.) This information should be sufficient to write an input deck for the
device simulator.
2. Every file created must start with the instance name. This provides a unique
way for identifying each file and avoids clobbering other files.
The npiscwritestanford.c file can be used as a starting point and modified for the needs of the
new device simulator. The code is well documented and should not be difficult to follow.
npiscreadMYSIMULATOR.c
This routine is responsible for reading the results of the device simulation. When creating the input
deck, one should have specified output files for saving the current vector and conductance matrix of
the device solution. Given the pointer to the bias data structure which contains the instance name for
identification and given a pointer to a solution data structure, this routine reads the simulation results
and stores them in the solution data structure. It must also multiply by the area factor if it is not done
during the device simulation.
The npiscreadstanford.c file can be used as a starting point and modified for the needs of the new
device simulator. The code is well documented and should not be difficult to follow.
npiscmethod.c
In PISCES, there is a card called “method” that is used to define the numerical methods during the
device simulation. This method is the same for all device simulations for a specific model. Hence, the
mixed-mode interface has been designed to accept this method card in two formats (Refer to “Method
and Model Parameters” on page 288 for more information.), store in its internal structures, and pass it
to the npiscwriteMYSIMULATOR.c routine through the bias data structure. Since other device
simulators have something similar, this routine is modified to take that into account.
There are two ways to determine the method line. The easiest way to implement the method card is to
read it from a file named [Model Name].method which is the default in npiscmethod.c.
In order to eliminate an extra file, a binary encoded method parameter is added to the SPICE .model
parameters. Essentially, a binary number (written in decimal format on the .model card) determines
whether a logical expression on the PISCES method line is true or false. The numdevdefs.h file
contains defined variables for all the possible actions. Each action has a number in base two that
corresponds to a bit in a long integer.
The npiscmethod.c routine takes the integer value given on the .model card and converts that to
character string for the numerical device input deck. This file and the numdevdefs.h file are well
documented and the programmer should not have a problem understanding how to modify this routine.
If there is no method line, everything is initialized to NULL.
npiscmodel.c
This file is analogous to the npiscmethod.c file except it does the processing for the model card in
the numerical device input deck. Refer to this routine and the numdevdefs.h file for more
information.
After all changes have been made, the new files can be added to the makedefs file so that they are
compiled and linked when SPICE is built. Once included, the new device simulator can be used in a
mixed-mode simulation simply by specifying the proper level on the .model card in the SPICE
netlist.
References
[1] Yu, Zhiping, Robert W. Dutton, and Hui Wang. “A Modularized, Mixed IC Device/ Circuit
Simulation System.” Proceedings of the Synthesis and Simulation Meeting and International
Exchange. Kobe, Japan. April 1992.
[3] Mayaram, Kartikeya and Donald O. Pederson. “Coupling Algorithms for Mixed-Level Circuit
and Device Simulation.” IEEE Transactions on Computer Aided Design. Vol. II, No 8. August
1992. pp 1003-1012.
[4] Rollins, Gregory J. and John Choma. “Mixed-Mode PISCES-SPICE Coupled Circuit and
Device Solver.” IEEE Transactions on Computer Aided Design. Vol. 7, No. 8, August 1988. pp
862-867.
[5] Geist, Al, Adam Beguelin, Jack Dongarra, Weicheng Jiang, Robert Manchek, and Vaidy
Sunderam. PVM 3 User’s Guide and Reference Manual. Oak Ridge: Oak Ridge National
Laboratory, 1993.
[6] Johnson, B., T. Quarles, A. R. Newton, D.O Pederson, and A. Sangiovanni-Vincentelli. SPICE3
Version 3f User’s Manual. Berkeley: Regents of the University of California, 1992.
[7] Dutton, R. W., F. Rotella, Z. Sahul, L. So, and Z. Yu. “Integrated TCAD for OEIC
Applications.” International Symposium on Optoelectronics for Information and Microwave
Systems. Proceedings of the SPIE - The International Society for Optical Engineering. Los
Angeles, CA. January 1994.
[8] Quarles, Thomas L. “Adding Devices to SPICE3.” Memorandum No. UCB/ERL M89/45.
Berkeley: University of California, 1989.
[9] Quarles, Thomas L. “SPICE3 Implementation Guide.” Memorandum No. UCB/ERL M89/44.
Berkeley: University of California, 1989.