蒙地卡羅解波茲曼方程式演算法準確性再提升以及平行模擬

(1)

行政院國家科學委員會補助專題研究計畫成果報告

※※※※※※※※※※※※※※※※※※※※※※※※※

※ ※

※ 蒙地卡羅解波茲曼方程式演算法準確性再提升 ※

※ 以及平行模擬 ※

※ ※

※※※※※※※※※※※※※※※※※※※※※※※※※

計畫類別：個別型計畫

計畫編號：NSC 90-2118-M-009-012

執行期間：90 年 8 月 1 日至 91 年 7 月 31 日

計畫主持人：周幼珍

本成果報告包括以下應繳交之附件：

□赴國外出差或研習心得報告一份

□赴大陸地區出差或研習心得報告一份

□出席國際學術會議心得報告及發表之論文各一份

□國際合作研究計畫國外研究報告書一份

執行單位：交通大學統計所

中華民國 91 年 10 月 31 日

(2)

蒙地卡羅解波茲曼方程式演算法準確性

再提升以及平行模擬

*

計畫編號：NSC 90-2118-M-009-012

執行期限：90 年 8 月 1 日至 91 年 7 月 31 日

主持人：周幼珍博士交通大學統計所

一、中文摘要波茲曼運動方程式(Boltzman transport equation, BTE)描述半導體微觀的狀況已被廣泛運用與研究。模擬電子運動過程中，蒙地卡羅法(Monte Carlo, MC)解耦合的泊松(Poisson)-BTE 模型不僅複雜、費時且有統計雜訊問題。本研究運用適應性網格切割法提升求解上的穩定性與效率性。測試結果發現此方法兼具收斂性及準確性，同時也可以有效實行與平行計算環境上。關鍵詞：波茲曼運動方程式、蒙地卡羅模擬法、模擬收歛性、平行計算 Abstract

We propose a parallel and adaptive 1-irregular mesh solution method for the MC solution of semiconductor Poisson-BTE transport model. We iteratively solve the coupled Poisson-BTE transport equations using MC particle method on an unstructured mesh for 1D/2D n+-n-n+ device. Most of conventional MC approaches applied to solve semiconductor Poisson-BTE transport model rely on fine structure mesh, so it is a time-consuming task and has significantly statistical noise. MC simulation on an unstructured 1-irregular mesh has more computational efficiency than that on structured mesh. It can be utilized for large-scale and 2D/3D nanodevice simulation. Computational results are presented to show the novelty and robustness of this method. Keywords: BTE, MC method, Convergence,

Parallelization

二、緣由與目的

Advances in semiconductor fabrication technology have allowed the continual reduction of minimum feature sizes consistent with Moore's Law. Gordon Moore made his famous observation in 1965, just 4 years after the first planar integrated circuit (IC) was discovered. The press called it "Moore's Law" and the name has stuck. In his paper [2], Moore predicted that the number of transistors per IC would double every 18 months. The Semiconductor Industry Association has projected this trend to continue for at least the next 15 years, aggressively pushing the bit count of DRAM by a factor of 4 every 3 years.

Among the grand challenges for high performance computing, nanoscale device physics, simulation, and characterization are of particular interest for their roles in developing a key industrial technology. This continual downscaling of device sizes has placed tremendous demands on predicting the behavior of designs [3,4]. TCAD/ECAD tools [5] for device simulation can considerably reduce the development cost and turnaround time for new devices. However, aggressive scaling of semiconductor technology is pushing the capabilities of today's simulation tools. Carrier transport in commercial TCAD packages is usually described with the drift diffusion or hydrodynamic models. However, these models may not adequate, because they hold only if the carriers are supposed to be in thermal equilibrium with the lattice (which is not the case in the high field and almost ballistic situations encountered in deep- submicron and nanoscale devices). The Monte Carlo Method (MCM) for semiconductor device simulation [6-12], which is based on a microscopic approach of

*This report is based on our paper [1] “Adaptive Monte Carlo Method for the Solution of Semiconductor Poisson-Boltzmann Transport Equations” in WSEAS

(3)

2

electrons in semiconductors, offers an alternative that is becoming increasingly important as device technology now focuses on the nanometer scale ULSI regime, where device sizes are less than 0.18 µm. From the MCM simulation, a detailed description of the underlying physics of the device behavior can be gleaned. Nevertheless, the large computational requirements of MCM, which may be on the order of hours or days of CPU time even for 2D treatments on high performance systems, prohibit its direct use in practical device engineering. An efficient application of MCM has therefore remained either for research or as a method of parameter extraction for higher-level models.

In this work, based our previous works for deep- submicron semiconductor device simulation with adaptive and parallelization methods, we propose an adaptive 1-irregular mesh simulation technique for MCM solution of the Poisson-Boltzmann Transport Equations (P-BTE) on a Linux-cluster. This approach is a direct extension of the parallel adaptive computing method for solving device models in our earlier works [13-19]. We for the first time iteratively solve the coupled P-BTE model using MCM on an unstructured mesh for 1D/2D n+-n-n+ device. Conventional MCMs for solving P-BTE model are with fine structure mesh. They are time-consuming task and have significantly statistical noise. Our adaptive MCM bases on adaptive 1-irregular mesh, monotone iterative, and a posteriori error estimation methods on a Linux-cluster. Our simulation results for an n+-n-n+ structure show that perform MCM simulation on an unstructured 1-irregular mesh has more computational efficiency than the conventional simulation on structured mesh. It can be utilized for multidimensional nanometer scale device simulation.

三、研究之波茲曼運動方程式

Due to the ongoing miniaturization in the design of devices mesoscopic models using kinetic transport equations for a distribution function describing the state of the electron gas become more and more

relevant: (1) the effects of impact-ionization, (2) tunneling, and (3) hot-carrier transport play an important role in the design of reliable small-scale devices. Therefore it is of paramount importance to describe these phenomena with a TCAD tool [5].

The semiclassical BTE with the Poisson equation for the electric field is the starting point for a rigorous theory. The analytical solution of such a system is impossible from the practical point of view, and for this reason the MCM numerical algorithm has been developed. The MCM is able to describe the behavior of small semiconductor devices even far from thermal equilibrium, because it does not make any restrictive artificial assumptions on the dynamics. Furthermore, the full band structure of the semiconductor and the physical scattering rates can be taken into account. Simulating nanometer size transistors requires solving an integro-differential 5D BTE coupled with the Poisson equation for electron transport. The BTE can be express as follows [10,12].

, ) ( coll t f k f F r f v t f δ δ = ∂ ∂ ⋅ + ∂ ∂ ⋅ + ∂ ∂ r h r r r ₍₁₎

where f = f(rr,kr,t) is the particle density distribution function in the momentum space at time t and location rr=(r_x,r_y,r_z). The

) , , (v_x v_y v_z

vr= is the particle velocity, rr is the position coordinate in phase space,

) , , (k_x k_y k_z

kr= is the crystal momentum

coordinates in phase space. The term F rr(r)

is the external force on the particle at location

rr. The collision term on the right hand side is ,' )} ,' ( ,.)] (., 1 ,.)[ ' (., ) ' , ( ,.)] ' (., 1 ,.)[ (., { ) ( k d k k S k f k -f k k S k f k f t f coll r r r r r r r r r ⋅ − ⋅ − − =

∫

δ δ (2)

where the S(kr,'kr) is the particle transition rate due to scattering mechanism. In many ways, it is similar to the BTE of neutron transport as well as rarefied gas dynamics, but for electrons transport in nanoscale structures, there are some added

(4)

complexities.

n+ _n+

n

x y

Fig. 1. A 2D cross section view for a n+-n-n+ structure

As shown in Fig. 1, the P-BTE model is subject to mixed type boundary conditions in 1D or 2D a simulation domain. On the source and drain contacts, the Dirichlet type boundary condition is applied for the electrostatic potential. We assume the boundary conditions for particles obeying the conservation law.

Both of the micro- and macro-scopic quantities, such as the distribution function, mean velocity, and temperature can be evaluated by statistical methods. The computed distribution function f can be used for the macroscopic physical quantities calculation. We apply the following formula for macroscopic electron density calculation

, ) , , ( ) , (r t =

∫

f r k t dk n r r r r (3)

use Eq. (4) to calculate the average electron velocity , ) , , ( 1 ) , ( >=

∫

< vf r k t dk n t r vr r r r r r (4) and follow Eq. (5) to calculate the average electron energy , ) , , ( 1 ) , ( >=

∫

< Ef r k t dk n t r E r r r r (5) where n is solved with Eq. (3).

In the simulation we consider here the parabolic band structure _E₍_k₎ 2_kr2_/₂_m*

h r

=

and velocity of the particle is evaluated with

*

/ m

kr

h . The scattering mechanisms applied

here consist of: (1) acoustic phonon scattering, (2) optical phonon scattering, (3) ionized impurity scattering, and (4) virtual scattering [6]. The BTE is coupled with the Poisson equation for the electric field

), (n p D q s∇ = − + ⋅ ∇ ε φ (6)

where the function D=−(N_D+ −N−_A) is the specified ionized net doping profile.

四、求解演算法

The conventional MCM consists of a simulation of the motion of electrons (and holes) inside the crystal, subject to the external and self-consistent electric field and of given scattering mechanisms [12]. The simulation starts with electrons in a given initial conditions for momentum and positions. During the free flight, the external forces are made to act according to the usual Newton's law. Then the scattering mechanism is chosen as responsible for the end of the free flight, according to the relative probabilities of all possible scattering mechanisms. From the differential cross section of this mechanism a new state after scattering is randomly chosen as initial state of the new free flight. The Poisson equation is solved to obtain the new electric field and the entire process is iteratively repeated. In this way a stochastic solution to the BTE is obtained. By running this procedure, the carrier history of each particle is recorded and one obtains information on some quantity, such as the distribution function, drift velocity, mean energy, and temperature stress tensor by taking the average of that quantity during the carrier histories in a time interval.

Various P-BTE device modeling and simulation have been developed for the approximation solution on a structured mesh and have their advantages [6-12]. However, the extensive computation required makes it

(5)

4

impractical for a device design on a regular basis. Concerning the computational efficiency, we present our adaptive MCM for the P-BTE model simulation. We state the adaptive computing algorithm and the parallelization procedure for the P-BTE Monte Carlo simulation. It relies on the estimations of the solution gradient and variation of carrier lateral current density along the device channel surface. We outline the simulation procedure for the MCM simulation on 1-irregular mesh.

The adaptive computing algorithm applied in Step 9 above is followed our earlier works [13-19] and is summarized in Alg. 2. We utilize this algorithm to solve the Poisson equation in 1D/2D structure.

If the estimator is less than a preset error tolerance, the process will be terminated and the approximated solution can be post-processed for next iterations. Otherwise, our refinement scheme will refine the mesh using the maximum gradient of electrostatic potential φ and the variation of current density Jn as error estimation. The proposed

parallel computing procedure for accelerating the MCM execution speed is outlined in the following algorithm. Alg. 3 is based on the partition of total particles numbers.

五、結果與討論

We present simulated results; the simulation example is a 1D 0.4µm n+-n-n+ semiconductor structure with 0.2 µm source and drain contacts [3]. The applied voltage VD is fixed at 2 V. The simulation consists of

50000 particles more than 100000 iteration loops to reach the self-consistent solution.

Fig. 2 shows the simulated electron density, where the solid line is the specified doping profile and the dot one is the corresponding electron density in log scale. Fig. 3 shows the electron density variation between the junctions. Our adaptive scheme successfully locates the solution variation. We plot the potential distribution and electric field in Figs. 4 and 5, respectively. Compared

Algorithm 1. (Adaptive single-particle Monte

Carlo simulation method for P-BTE model)

Step 1. Compute the maximum of total scattering

rate.

Step 2. Initialize the physical quantities for all

particles.

Step 3. Determine the free flight time of the

particle.

Step 4. Determine each particle location and its

electric field with unstructured mesh data.

Step 5.Calculate each particle’s momentum and

energy before the collision.

Step 6. Determine the scattering process.

Step 7. Compute motion angles θ and φ after the

collision.

Step 8. Calculate their momentum and energy after

the collision.

Step 9. After performed the simulation for all

particles, we solve the Poisson equation with the adaptive computing algorithm. Otherwise return to Step 3.

Step 10. Update the electric field for all

unstructured mesh.

Step 11. Stop the iteration loops and perform the

post- processing, if the number of iterations is grater than a specified stopping iteration counts.

Algorithm 2. (Adaptive computing algorithm for

the Poisson equation.)

Step 1. Poisson Model formulation Step 2. Domain discretization. Step 3. Finite volume approximation.

Step 4. System of the corresponding nonlinear

algebraic equations: Ax = -F(x).

Step 5. The MI nonlinear solver. Step 6. A posteriori error estimation. Step 7. Mesh refinement.

Algorithm 3. (Parallel computational procedure)

Step 1. Initialize the MPI and configurations. Step 2. Setup unstructured mesh dynamical array

structure.

Step 3. Count number of particles and apply a dynamic MCM partition algorithm to

determinate number of processors in the simulation.

Algorithm 3.1. (dynamic MCM partition

algorithm)

Step 1. Count the number of total particles. Step 2. Find out the optimal number of

processors.

Step 3. Determine how many particles should

be assigned to each processor.

Step 4. Along x- or y-direction in device

domain, search (from left to right and bottom to top) and assign particles to these processors sequentially.

Step 4. All assigned CPUs perform the MCM

(6)

with the Fig. 10, the computed electric field provides a very good error control scheme and mesh refinement indicator.

Distance (m)

0 2e-7 4e-7 6e-7 8e-7

E lec tro n den si ty ( m ^3 ) 1e+21 1e+22 1e+23 1e+24 1e+25 Doping profile Electron density

Fig. 2. The device doping profile and its electron density at VD = 2.0 V.

Distance (m)