Semiconductor process device simulation method and storage medium storing simulation program

[0002] The present invention relates to a semiconductor process device simulation method of simulating the manufacturing process and electrical characteristics of a semiconductor device using a computer, and a storage medium storing a simulation program and, more particularly, to a simulation method of efficiently and numerically solving, on a computer, simultaneous multi-dimensional linear equations based on physical laws used for a simulation.

[0003] As semiconductor integrated circuits such as a VLSI have become smaller, manufacturing has become increasingly complicated. As a result low-temperature processing, precise annealing processing, and accurate process design are of increasing importance.

[0004] On such a background, in recent years, device simulations are generally used for development and research of, e.g., the manufacture of semiconductor devices. A factor for this is a great progress in computers which actually execute the device simulation. Additionally, more convenient interfaces for the computers are becoming available. The progress in the analysis technique allows simulations for objects in a wider range and even analysis of breakdown phenomena or current concentration in a semiconductor device.

[0005] In computer simulations of thermal impurity diffusion which is one of the manufacturing processes of a semiconductor device described in reference 1: Ryo Dan, “Process Device Simulation Technology”, pp. 26-28, first a region to be analyzed is divided into meshes, and a diffusion equation is discretely defined for each mesh point. The diffusion equations are converted into linear equations and then into simultaneous linear equations by Newton's method or the like so as to solve the diffusion equations.

[0006] The electrical characteristics of a semiconductor device prepared by manufacturing processes including thermal impurity diffusion, ion implantation, and thermal oxidation can be simulated by the method described in reference 2 (Dan Ryo, “Process Device Simulation Technology”, pp. 105-134), according to which a region to be analyzed is divided into meshes, and a Poisson equation and current continuous equation are discretely defined for each mesh point. These equations are converted into linear equations and then into simultaneous linear equations by Newton's method or the like, and the equations are solved on a computer, thereby simulating the electrical characteristics of the semiconductor device.

[0007] In these semiconductor process device simulation methods, it is important in practice to solve large scale simultaneous linear equations having as many dimensions as the number of mesh points using a computer with as small a memory and as high a speed as possible.

[0008] The function of solving these simultaneous linear equations is called a “matrix solver” for which various methods have been proposed. For example, for an equation having four dimensions factorization may be used.

[0009] And, for an equation having a large number of dimensions, an iterative method is used because of limitations on the memory capacity and calculation time needed. When the coefficient matrix of simultaneous linear equations to be solved is symmetrical, the ICCG (Incomplete Choleske and Conjugate Gradient) method is used. This is a CG (Conjugate-Gradient) method with preprocessing.

[0010] When the coefficient matrix is asymmetrical, a series of techniques called the Krylov subspace method are generally used because of advantages in convergence.

[0011] With both the ICCG and Krylov subspace methods, forming an algorithm for multiplying the original coefficient matrix with an approximate inverse matrix of the coefficient matrix of the equation to be solved permits a reduction in the number of conditions of the coefficient matrix, whereby the approximate inverse matrix of the coefficient matrix of the equation may be solved by a reduced number of iterations. This operation is called “preconditioning”, as described above. The approximate inverse matrix of the coefficient matrix is called a “preconditioning matrix”.

[0012] As the inverse matrix of the preconditioning matrix, a matrix obtained by “incomplete LU-factorization”, i.e., LU-factorization, with limited fill-in generation positions is often used. Fill-in is a process of LU-factorization represented by equations (1) and (2) below and means that at a position (i, j) where an element value A_ijis 0 in the original coefficient matrix, a new element L_ijor U_ijother than 0 is generated. When the generated fill-in is not rejected, and taken into consideration in the subsequent process of $\begin{array}{cc}{L}_{ij}=A_{ij}-\sum _{k=1}^{j-1}{L}_{ik}{U}_{kj}\left(k<j<i\right)& \left(1\right)\\ U_{ij}=L_{ii}^{-1}\left(A_{ij}-\sum _{k=1}^{i-1}{L}_{ik}{U}_{kj}\right)\left(k<i<j\right)& \left(2\right)\end{array}$

[0013] LU-factorization, conventional complete LU-factorization is performed.

[0014] In preprocessing based on incomplete LU-factorization, it is required to perform optimum processing which satisfies the following two contradictory conditions. First, fill-in generation should be minimized to suppress the amount of memory used and calculation time needed for performing the calculation. Second, a result as close as possible to complete LU-factorization should be obtained to increase convergence of the iterative method.

[0015] A conventional method of efficiently performing such processing is described in reference 3 (Shin Nakamura and Akio Nakagawa, High-speed Iterative Method for Two-dimensional Bipolar MOS Composite Device Simulator TONADDEII”, IEICE Technical Report, pp. 64-65) or reference 4 (Megumi Kawakami, Masahiro Sugaya, and Shiroo Kamohara, “A New High-speed Non-equilibrium Point Defect Model for Annealing Simulation”, SISPAD '96, p. 94, FIG. 1a) is used.

[0016] In the method described in reference 3 or 4, as shown in the flow chart of FIG. 8, first, the maximum value n of the number of equations to be solved which are defined for each point is obtained in step 801. In step 802, the rows and columns of the coefficient matrix are put into a group in units of square submatrices of n×n for each mesh point, thereby forming blocks of coefficient matrices. In step 803, while virtually regarding an n×n square submatrix as one matrix, incomplete LU-factorization processing based on calculation of the square submatrices is performed. FIG. 9 shows the structure of the coefficient matrix formed by the above operation.

[0017] With the operation in units of square submatrices, incomplete LU-factorization close to complete LU-factorization can be more effectively performed using the same memory capacity as would be in processing units of scalar elements. This is because the fill-in generation pattern is determined while giving prominence to a combination of equations at one mesh point.

[0018] Objects to which the above-described conventional semiconductor process device simulation technology is applied are disclosed in Japanese Patent Laid-Open No. 6-53155 or 8-213334.

[0019] However, the above-described conventional semiconductor process device simulation technology has the following problem. That is, when a equation is derived from a certain physical law holds for only part of an analysis region to be processed, or no equation is defined, excess memory capacity is required to process incomplete LU-factorization.

[0020] This is because in making the computer calculations, a square submatrix having a uniform size of n×n is used as the processing unit of incomplete LU-factorization. For example, in the above-described simulation, the submatrix size of n×n is maintained even at a mesh point where n equations are not defined. The computer performing the calculations must form an overall coefficient matrix while inserting “1” to the corresponding principal diagonal portions of the submatrix and set the coefficient matrix on the memory. For this reason, excess memory capacity is used.

[0021] When the computer executes its calculation, excess memory is accessed. This wasteful processing adds to processing time.

[0022] It is the main object of the present invention to provide high speed simulation without using an excess memory for calculation.

[0023] It is another object of the present invention to provide a semiconductor process device simulation method which allows high-speed incomplete LU-factorization without using excess memory capacity even when an equation derived from a certain physical law and used to simulate the manufacturing process or electrical characteristics of a semiconductor device using a computer holds at only one portion of a region to be analyzed, or no equation is defined, and a storage medium storing a simulation program therefor.

[0024] In order to achieve the above objects, according to one aspect of the present invention there is provided a semiconductor device manufacturing process simulation method to aid manufacturers in forecasting electrical characteristics of semiconductor devices by performing a plurality of matrix manipulations of terms representing physical properties of the semiconductor devices, the matrices to be manipulated representing multidimensional simultaneous linear equations that are to be solved by a matrix solver that uses an iterative method in which the matrices are preconditioned by incomplete LU-factorization, the method for generating fill-ins for the matrices comprising the steps of: a first step of dividing the surface of a semiconductor device to be processed into a plurality of rectangles forming a mesh of predetermined size; a second step of assigning a numerical value to each mesh point of the mesh; a third step of setting equations representing a relationship among the plurality of numerical values; a fourth step of generating a coefficient matrix constituted by a plurality of principal diagonal submatrices each of which is arranged at each one of principal diagonal positions corresponding to each mesh point and representing a self feedback function to the mesh point, the coefficient matrix having rows and columns in numbers corresponding to the number of mesh points, and a plurality of non-principal diagonal submatrices each of which is arranged on any one of rows and columns passing through the principal diagonal positions corresponding to the mesh point and representing an interaction between the mesh point corresponding to the principal diagonal located on the same row or on the same column of the coefficient matrix with the non-principal diagonal submatrix and an adjacent mesh point connected to the mesh point through a mesh branch; a fifth step of performing calculation for the submatrices while regarding each submatrix of the coefficient matrix as one element, thereby performing incomplete LU-factorization of the coefficient matrix, wherein each of the principal diagonal submatrices is a square having rows and columns equal in number to equations set for a mesh point corresponding to the principal diagonal submatrix, each of the non-principal diagonal submatrices arranged in a row direction in correspondence with each of the mesh points is a matrix having rows equal in number to equations set at a mesh point corresponding to the principal diagonal submatrix located in the row and columns equal in number to equations set at an adjacent mesh point connected to the mesh point through a mesh branch, and each of the non-principal diagonal submatrices arranged in a column direction in correspondence with each of the mesh points is a matrix having columns equal in number to equations set at a mesh point corresponding to the principal diagonal submatrix located in the column and rows equal in number to the equations set at an adjacent mesh point connected to the mesh point through a mesh branch; and a sixth step of producing a signal indicative of the result of said calculation.

[0025] In another aspect of the invention there is provided a computer readable memory storing a semiconductor device manufacturing process simulation program to aid manufacturers in forecasting electrical characteristics of semiconductor devices by performing a plurality of matrix manipulations of terms representing physical properties of the semiconductor devices, the matrices to be manipulated representing multidimensional simultaneous linear equations that are to be solved by a matrix solver that uses an iterative method in which the matrices are preconditioned by incomplete LU-factorization, the program including a routine for generating fill-ins for the matrices, the routine comprising the steps of: causing a computer to perform the following functions: a first function of dividing a surface of a semiconductor device to be processed into a mesh of predetermined size; a second function of assigning a numerical value to each mesh point of the mesh; a third function of setting equations representing a relationship among the plurality of numerical values; a fourth function of generating a coefficient matrix constituted by a plurality of principal diagonal submatrices each of which is arranged at each one of principal diagonal positions corresponding to each mesh point and representing a self feedback function at the mesh point, the coefficient matrix having rows and columns in numbers corresponding to the number of mesh points, and a plurality of non-principal diagonal submatrices each of which is arranged on any one of rows and columns and representing an interaction between the mesh point corresponding to the principal diagonal positions located on the same row or on the same column of the coefficient matrix with the non-principle diagonal submatrix and an adjacent mesh point connected to the mesh point through a mesh branch; and a fifth function of performing calculation for the submatrices while regarding each submatrix of the coefficient matrix as one element; thereby performing incomplete LU-factorization of the coefficient matrix, wherein each of the principal diagonal submatrices is a square matrix having rows and columns equal in number to equations set forth for a mesh point corresponding to the principal diagonal submatrix, each of the non-principal diagonal submatrices arranged in a row direction in correspondence with each of the mesh points is a matrix having rows equal in number to equations set at a mesh point corresponding to the principal diagonal submatrix and located in the row and columns equal in number to equations set at an adjacent mesh point connected to the mesh point through a mesh branch, and each of the non-principal diagonal submatrices arranged in a column direction in correspondence with each of the mesh points is a matrix having columns equal in number to equations set a mesh point corresponding tot he principal diagonal submatrix located in the column and rows equal in number to equations set at an adjacent mesh point connected to the mesh point through a mesh branch.

[0026] With this arrangement, the non-principal diagonal submatrix not always will be a square matrix. When the number of equations set at an adjacent mesh point connected to a mesh point of interest through a mesh branch is smaller than that set for the mesh point of interest, a non-square matrix is obtained, so the number of calculations is reduced as compared to a case in which a square matrix is simply assigned.

[0027] FIG. 1 is a block diagram showing the arrangement of a simulation apparatus for realizing a semiconductor process device simulation method according to the present invention;

[0028] FIG. 2 is a flow chart showing the flow of basic operation of semiconductor process device simulation executed by the simulation apparatus (CPU 100) shown in FIG. 1;

[0029] FIG. 3 is an explanatory view showing an example of the structure of a coefficient matrix formed in the simulation operation;

[0030] FIG. 4 is an explanatory view schematically showing the section of an object to be simulated in the first embodiment of the present invention;

[0031] FIG. 5 is an explanatory view showing a coefficient matrix used to simulate the object to be simulated, shown in FIG. 4;

[0032] FIG. 6 is an explanatory view schematically showing the section of an object to be simulated in the second embodiment of the present invention;

[0033] FIG. 7 is an explanatory view showing a coefficient matrix used to simulate the object to be simulated, shown in FIG. 6;

[0034] FIG. 8 is a flow chart showing the flow of basic operation of conventional semiconductor process device simulation; and

[0035] FIG. 9 is an explanatory view showing an example of the structure of a coefficient matrix formed in a conventional simulation operation.

[0036] The present invention will be described below with reference to the accompanying drawings.

[0037] FIG. 1 shows the arrangement of a simulation apparatus for realizing a semiconductor process device simulation method according to the present invention.

[0038] This simulation apparatus comprises a central processing unit (CPU) 100, a memory 102, a display unit 103, an input unit 104, and an external storage means 105 which are connected through a bus 106. More specifically, the simulation apparatus is realized by a workstation, a personal computer, or any other computer system. For example, as the external storage means 105, a magnetic disk, a semiconductor memory, or any other general storage medium can be used.

[0039] In this simulation apparatus, a program (simulation program) for simulation is stored in the external storage means 105. When the operator inputs an instruction for the start of simulation from the input unit 104, the CPU 101 loads the program from the external storage means 105 into a predetermined area in the memory 102. The operator prepares data necessary for simulation (matrix) in the memory 102 by, e.g., inputting the data from the input unit 104. The data necessary for simulation may be stored in the external storage means 105 in advance and loaded in the memory 102. The CPU 101 performs predetermined calculations using the data in accordance with the program loaded in the memory 102, thereby performing predetermined simulation. Finally, the CPU 101 displays the simulation result on the display unit 103.

[0040] FIG. 2 shows the flow of basic operation of semiconductor process device simulation executed by the simulation apparatus (CPU 101) shown in FIG. 1.

[0041] The operation will be briefly described. Beginning with steps 201 to 203 first, in step 201, a mesh having a predetermined pitch is set on the object to be simulated by, e.g., inputting data from the input unit 104 (FIG. 1), and stored in the memory 102. Assume that I mesh points are set. In step 202, a plurality of values are defined for each mesh point by, e.g., inputting them from the input unit 104 and stored in the memory 102. In step 203 equations representing the relationship between the plurality of values are formed and stored in the memory 102.

[0042] In step 204, the CPU 101 checks the number n(i) of equations to be solved, which represent the relationship between the plurality of values defined for each mesh point i (i is an integer, 1≦i≦I) and are stored in the memory 102, on the basis of the loaded program, and stores the value in the memory 102 shown in FIG. 1.

[0043] In step 205, the CPU 101 puts equations to be solved, which are stored in the memory 102, into groups in units of mesh points, and forms, in the memory 102, a coefficient matrix represented as blocks arranged at nonuniform pitches and having elements at the ith row and jth column as an n(i)×n(j) submatrix. More specifically, a coefficient matrix having rows and columns equal in number to mesh points is set. A principal diagonal submatrix representing the self feedback function at each mesh point is arranged at a portion (element) corresponding to the mesh point. In addition, at one of the row and column passing through the principal diagonal position corresponding to each of mesh points other than the principal diagonal of the coefficient matrix, a non-principal diagonal submatrix representing the interaction between a mesh point corresponding to the principal diagonal position and an adjacent mesh point connected to the mesh point through a mesh branch is arranged. The principal diagonal submatrix is a square matrix having rows and columns, the number of which is equal to the number n(i) of equations defined for the corresponding mesh point. A non-principal diagonal submatrix arranged in the row direction in correspondence with each mesh point is a matrix having rows the number of which is equal to the number n(i) of equations defined for the mesh point and columns the number of which is equal to the number n(j) of equations defined for an adjacent mesh point connected to the mesh point through a mesh branch. A non-principal diagonal submatrix arranged in the column direction in correspondence with each mesh point is a matrix having columns the number of which is equal to the number n(i) of equations defined for the mesh point and rows the number of which is equal to the number n(j) of equations defined for an adjacent mesh point connected to the mesh point through a mesh branch.

[0044] Finally, in step 206, the CPU 101 performs calculations for the submatrices while regarding the submatrix at the ith row and jth column, which is formed on the memory 102, as one element, thereby performing incomplete LU-factorization of the coefficient matrix.

[0045] The simulation operation of the present invention after the equations representing the relationship between the values assigned for each mesh point are set will be described next with reference to the flow chart of FIG. 2 and the view of FIG. 3. FIG. 3 shows an example of the structure of the coefficient matrix formed in the simulation operation.

[0046] The number of equations at each mesh point on the surface of the object to be simulated, i.e., a semiconductor device on which a mesh is set, is checked (FIG. 2: step 204). In the following description, assume that a total of seven mesh points are present, and three equations are defined at mesh point 1, two equations are defined at mesh point 2, one equation is defined at mesh point 3, three equations are defined at mesh point 4, two equations are defined at mesh point 5, one equation is defined at mesh point 6, and four equations are defined at mesh point 7 so as to represent the states at the respective mesh points. That is, n(1)=3, n(2)=2, n(3)=1, n(4)=3, n(5)=2, n(6)=1, and n(7)=4.

[0047] The equations to be solved are put into groups in units of mesh points to form coefficient matrix blocks (FIG. 2: step 205). When the above-described mesh is set, coefficient matrices as shown in FIG. 3 are obtained. Referring to FIG. 3, a hatched portion enclosed by bold solid lines represents a block, and broken lines represent boundaries between scalar elements in each block. For example, the block at the seventh row and fourth column forms a 4×3 non-square matrix.

[0048] A submatrix is regarded as one element. The equations are solved by performing incomplete LU-factorization according to the above-described equations (1) and (2) by calculation of the square and non-square matrices, thereby simulating the state at each mesh point (step 206). For example, in FIG. 3, an element L₃₃is generated by factorizing the third row and third column by subtracting the product of a 1×3 submatrix L₃₁and a 3×1 submatrix U₁₃, which have been decomposed, from a 1×1 coefficient matrix element A₃₃, i.e., “L₃₃=A₃₃−L₃₁U₁₃. . . (3).”

[0049] The present invention will be described using a more specific example. An example in which the diffusion state of an impurity doped into a silicon substrate having an oxide film is simulated will be described.

[0050] In this impurity diffusion simulation, a mesh having a predetermined pitch is allocated to a section (two-dimensional region) of a silicon substrate having an oxide film, and predetermined equations are set for each mesh point. FIG. 4 shows the two-dimensional region. A region 401 divided into eight blocks is the region of the silicon substrate. A region 402 divided into two blocks on the region 401 indicates the oxide film formed on the silicon substrate. A line segment 403 represents the boundary between the silicon substrate and the oxide film.

[0051] Mesh points 404 having numbers “10” to “21” and marked with circles in the region 401 are mesh points assigned in the two-dimensional region of the silicon substrate. Mesh points 405 having numbers “1” to “3” and marked with triangles in the region 402 are mesh points assigned on the surface of the two-dimensional oxide film region. Mesh points 406 having numbers “7” to “9” and marked with rectangles are mesh points representing the silicon substrate side of the interface between the silicon substrate and the oxide film. Mesh points 407 having numbers “4” to “6” and marked with stars are mesh points representing the oxide film side of the interface between the silicon substrate and the oxide film. The mesh points are coupled through mesh branches 408.

[0052] In this impurity diffusion simulation, diffusion simulation of boron is performed using a diffusion model paired with a point defect. In the model for simulation, the following physical phenomenon is used on the silicon substrate. First, assume that interstitial silicon and silicon vacancies as point defects are present in the silicon substrate, and the interstitial silicon and silicon vacancies diffuse while being recombined with each other. Boron as an impurity diffuses while forming pairs with interstitial silicon. When the pairs are dissociated, the boron atoms are set at the substitution positions of the crystal lattice of the silicon substrate.

[0053] According to the above point defect-paired diffusion model, five equations (4) to (8) below are equations to be solved at each of the mesh points 404 from “10” to “21” in the silicon substrate. $\begin{array}{cc}\begin{array}{c}\frac{\partial C_{BI}}{\partial t}=\nabla \left[D_{BI}\nabla C_{BI}+\frac{q{D}_{BI}}{kT}{Z}_{BI}{C}_{BI}\nabla \phi \right]+\\ k_{BI-f}{C}_{Ba}{C}_I-k_{BI-r}{C}_{BI}\end{array}& \left(4\right)\\ \frac{\partial C_{Ba}}{\partial t}=k_{BI-r}{C}_{BI}+k_{BI-f}{C}_{Ba}{C}_I& \left(5\right)\\ \begin{array}{c}\frac{\partial C_I}{\partial t}=\nabla \left[D_I\nabla C_I+\frac{q{D}_I}{kT}{Z}_I{C}_I\nabla \phi \right]+\\ k_{BI-r}{C_B}_I-k_{BI-f}{C}_{Ba}{C}_I-k_{IV}\left(C_I{C}_V-C_I^{*}{C}_V^{*}\right)\end{array}& \left(6\right)\\ \frac{\partial C_V}{\partial t}=\nabla \left[D_V\nabla C_V+\frac{q{D}_V}{kT}{Z}_V{C}_V\nabla \phi \right]-k_{IV}\left(C_I{C}_V-C_I^{*}{C}_V^{*}\right)& \left(7\right)\\ \begin{array}{c}\nabla \left(ϵ\nabla \phi \right)=-q[Z_{BI}{C}_{BI}-C_{Ba}+Z_I{C}_I+Z_V{C}_V-\\ n_i\exp \left(\frac{q\phi }{kT}\right)+n_i\exp \left(-\frac{q\phi }{kT}\right)]\end{array}& \left(8\right)\end{array}$

[0054] Of the above equations, equation (4) represents diffusion, generation, and dissociation of boron-interstitial silicon pairs. Equation (5) represents bond and dissociation of boron. Equation (6) represents diffusion, bond, dissociation, and disappearance of interstitial silicon. Equation (7) represents diffusion and disappearance of silicon vacancies. Equation (8) is a Poisson equation representing the relationship among charged point defects, boron concentration, and electrostatic potential.

[0055] In equations (4) to (6), C_BIis the volume concentration of boron-interstitial silicon pairs, C_Bais the volume concentration of active boron at the substitution positions, C_Iis the volume concentration of interstitial silicon, and C_Vis the volume concentration of silicon vacancies.

[0056] Z_BIis the number representing the charged state of boron-interstitial silicon pairs, Z_iis the number representing the charged state of interstitial silicon, and Z_Vis the number representing the charged state of silicon vacancies.

[0057] D_BIis the diffusion constant of boron-interstitial silicon pairs, D_Iis the diffusion constant of interstitial silicon, and D_Vis the diffusion constant of silicon vacancies.

[0058] In addition, k_BI−fis the formation constant of boron-interstitial silicon pairs, k_BI−ris the dissociation reaction rate constant, and k_IVis the recombination ratio of interstitial silicon and silicon vacancies.

[0059] Furthermore, q is the unit charge, ∈ is the dielectric constant of silicon, k is the Boltzman constant, T is the absolute temperature, φ is the electrostatic potential, and n_iis the intrinsic carrier density of silicon.

[0060] C*_Iis the volume density of interstitial silicon in the thermal equilibrium state, and C*_Vis the volume density of silicon vacancies in the thermal equilibrium state.

[0061] Diffusion of boron in the oxide film can be described using a normal model when boron itself directly diffuses without forming pairs. In addition, in the oxide film, the influence of the electrostatic potential can be disregarded.

[0062] Therefore, at the mesh points 405 in the oxide film, the concentration of boron as an impurity at these mesh points can be simulated by solving a normal diffusion equation such as equation (9) below: $\begin{array}{cc}\frac{\partial C_{Bc}}{\partial t}=\nabla \left[D_B\nabla C_{Bc}\right]& \left(9\right)\end{array}$

[0063] where C_Bcis the total concentration of boron in the oxide film, and D_Bis the diffusion constant of boron in the oxide film.

[0064] At the interface between the silicon substrate and the oxide film, a flux Js due to segregation of boron, which is given by equation (10) below, is present between the mesh point 406 and the mesh point 407: $\begin{array}{cc}\mathrm{Js}=h\left(C_{Bc}-\frac{C_{BI}}{m}\right)& \left(10\right)\end{array}$

[0065] where C_Bcis the total concentration of boron at the mesh point 407 on the oxide film side, C_BIis the concentration of boron-interstitial silicon pairs at the mesh point 406 on the silicon side, h is the transportation coefficient of boron at the interface between the silicon substrate and the oxide film, and m is the segregation coefficient of boron at the interface between the silicon substrate and the oxide film.

[0066] The interstitial silicon and silicon vacancies in the silicon substrate disappear at the interface between the silicon substrate and the oxide film. In consideration of this, equations to be solved at the mesh point 406 on the silicon substrate side in diffusion simulation of boron as an impurity are five equations (11) to (15) below: $\begin{array}{cc}\begin{array}{c}\frac{\partial C_{BI}}{\partial t}=\nabla \left[D_{BI}\nabla C_{BI}+\frac{q{D}_{BI}}{kT}{Z}_{BI}{C}_{BI}\nabla \phi \right]+\\ k_{BI-f}{C}_{Ba}{C}_I-k_{BI-r}{C}_{BI}+h\left(C_{Bc}-\frac{C_{Bt}}{m}\right)\end{array}& \left(11\right)\\ \frac{\partial C_{Ba}}{\partial t}=k_{BI-r}{C}_{BI}+k_{BI-f}{C}_{Ba}{C}_I& \left(12\right)\\ \begin{array}{c}\frac{\partial C_I}{\partial t}=\nabla \left[D_I\nabla C_I+\frac{q{D}_I}{kT}{Z}_I{C}_I\nabla \phi \right]+\\ k_{BI-r}{C_B}_I-k_{BI-f}{C}_{Ba}{C}_I-k_{IV}\left(C_I{C}_V-C_I^{*}{C}_V^{*}\right)-\\ k_{IS}\left(C_I-C_I^{*}\right)\end{array}& \left(13\right)\\ \begin{array}{c}\frac{\partial C_V}{\partial t}=\nabla \left[D_V\nabla C_V+\frac{q{D}_V}{kT}{Z}_V{C}_V\nabla \phi \right]-\\ k_{IV}\left(C_I{C}_V-C_I^{*}{C}_V^{*}\right)-k_{VS}\left(C_I-C_I^{*}\right)\end{array}& \left(14\right)\\ \begin{array}{c}\nabla \left(ϵ\nabla \phi \right)=-q[Z_{BI}{C}_{BI}-C_{Ba}+Z_I{C}_I+Z_V{C}_V-\\ n_i\exp \left(\frac{q\phi }{kT}\right)+n_i\exp \left(-\frac{q\phi }{kT}\right)]\end{array}& \left(15\right)\end{array}$

[0067] where k_ISof equation (13) is the recombination rate of the interstitial silicon at the interface between the silicon substrate and the oxide film, and k_VSin equation (14) is the recombination rate of silicon vacancies at the interface between the silicon substrate and the oxide film.

[0068] In addition, an equation to be solved at the mesh point 407 on the oxide film side is equation (16) below: $\begin{array}{cc}\frac{\partial C_{Bc}}{\partial t}=\nabla \left[D_B\nabla C_{Bc}\right]-h\left(C_{Bc}-\frac{C_{BI}}{m}\right)& \left(16\right)\end{array}$

[0069] The simulation operation of the first embodiment will be described with reference to the flow shown in FIG. 2. First, the number of equations defined at each mesh point is checked (step 204). In the first embodiment, at a mesh point such as the mesh point 404 or 406 on the silicon substrate side in FIG. 4, five equations such as equations (4) to (8) or (11) to (15) are defined. At a mesh point such as the mesh point 405 or 407 on the oxide film side, one equation such as equation (9) or (16) is defined. Subsequently, the equations to be solved are put into a group in units of mesh points to form a coefficient matrix represented by blocks shown in FIG. 5 (step 205). In FIG. 5, a coefficient matrix having 21 rows and 21 columns which correspond to numbers “1” to “21” assigned to the mesh points in FIG. 4 is obtained. In one row of the coefficient matrix shown in FIG. 5, all of the square submatrix blocks at the mesh points constituting the principal diagonals having the row numbers as the mesh point numbers shown in FIG. 4 and the square or non-square submatrix blocks at adjacent mesh points connected to the mesh point through the mesh branches 408 are arranged.

[0070] In the first embodiment, for a mesh point such as the mesh point 404 or 406 on the silicon side, the width of each block is 5, as indicated by row “7” and column “7” to row “21” and column “21”. To the contrary, for a mesh point such as the mesh point 405 or 407 to which numbers “1” to “6” are assigned on the oxide film side, the width of each block is 1, as indicated by row 1.111 and column “1” to row “6” and column “6”. In each block shown in FIG. 5, a blank portion represents that a numerical value to be calculated is not stored.

[0071] For example, for row “8” in FIG. 5, which corresponds to the mesh point 406 represented by “8” in FIG. 4, a 5×5 submatrix principal diagonal block of the mesh point 406 with “8” is present at row “8” and column “8”. This principal diagonal block has 20 numerical values to be calculated. As shown in FIG. 4, the mesh point 404 represented by “11” is present under the mesh point 406 represented by “8”, so a 5×5 submatrix block corresponding to the mesh point 404 represented by “11” and marked with a circle is present at row “8” and column “11” in FIG. 5. This submatrix block has seven numerical values to be calculated. As shown in FIG. 4, the mesh point 406 represented by “7” and the mesh point 406 represented by “9” are present on the left and right side of the mesh point 406 represented by “8”, respectively, so 5×5 submatrix blocks corresponding to these mesh points are present at row “8” and column “7” and row “8” and column “9”. Each of these submatrix blocks also has seven numerical values to be calculated. As shown in FIG. 4, the mesh point 407 represented by “5” is present above the mesh point 406 represented by “8”, so a 5×1 non-square submatrix block corresponding to the mesh point 407 is present at row “8” and column “5”, as shown in FIG. 5.

[0072] Similarly, for row “5” corresponding to the mesh point 407 represented by “5” in FIG. 4, a 1×1 submatrix principal diagonal block corresponding to the mesh point 407 is present at row “5” and column “5”. In correspondence with the mesh point 405 represented by “2” above the mesh point 407 represented by “5”, a 1×1 submatrix block corresponding to the mesh point 405 is present at row “5” and column “2” in FIG. 5. Since the mesh points 407 represented by “4” and “6” are present on the left and right sides of the mesh point 407 represented by “5 ”, 1×1submatrix blocks corresponding to the mesh points 407 are present at row “5” and column “14” in FIG. 5 and row “5” and column “6” in FIG. 5. In addition, since the mesh point 406 represented by “8” is present under the mesh point 407 represented by “5”, a 1×5 non-square submatrix block is present at row “5” and column “8” in FIG. 5.

[0073] Finally, each submatrix is regarded as one element, and calculation for square matrices and calculation for non-square matrices are performed to perform incomplete LU-factorization in accordance with equations (1) and (2), thereby solving the equations and simulating the states at the above-described mesh points (step 206).

[0074] As described above, in the first embodiment, calculation is performed using matrices whose calculation amount is smaller than that for a 5×5 submatrices. For this reason, wasteful calculation can be omitted, and the calculation time can be shortened. The principal diagonal block at each mesh point on the oxide film side is a 1×1 square matrix. A non-diagonal matrix at each mesh point at the interface is a 1×5 or 5×1 non-square matrix. Conventionally, since a 5×5 submatrix is formed even for these mesh points, an excess memory capacity is used, and wasteful calculation is performed.

[0075] The second embodiment of the present invention will be described next in detail. An example in which the electrical characteristics of a silicon substrate having an oxide film on the surface is simulated will be described below. For example, the MOS capacitor of a gate insulating film is simulated by calculating the potential value and carrier concentration at each lattice point of a mesh.

[0076] In this device simulation as well, a mesh having a predetermined pitch is assigned to a section (two-dimensional region) of the silicon substrate having the oxide film, and predetermined equations are set for each mesh point. FIG. 6 shows the two-dimensional region. A region 601 divided into eight blocks is the region of the silicon substrate. A region 602 divided into two blocks on the region 601 indicates the oxide film formed on the silicon substrate. A line segment 603 represents the boundary between the silicon substrate and the oxide film.

[0077] Mesh points 604 having numbers “4” to “18” and marked with circles in the region 601 are mesh points assigned in the two-dimensional region of the silicon substrate. Mesh points 605 having numbers “1” to “3” and marked with triangles in the region 602 are mesh points assigned on the surface of the two-dimensional oxide film region. Unlike the above-described impurity diffusion process simulation, only one mesh point having the nature on the silicon substrate side is assigned at the interface between the silicon substrate and the oxide film.

[0078] At each mesh point, equation (17) as a Poisson equation holds:

[0079] ∇(∈∇φ)=−q(N_D−N_A−n+p)  (17)

[0080] where ∈ is the dielectric constant of a substance, φ is the electrostatic potential, q is the unit charge, N_Dis the donor density, N_Ais the acceptor density, n is the electron density, and p is the hole density.

[0081] At a mesh point such as the mesh point 604 on the silicon substrate side, equations (18) and (19) as current continuous equations hold: $\begin{array}{cc}\frac{\partial n}{\partial t}=\nabla \left(D_n\nabla n-{\mu }_n\nabla \phi \right)+G-R& \left(18\right)\\ \frac{\partial p}{\partial t}=\nabla \left(D_p\nabla p-{\mu }_p\nabla \phi \right)+G-R& \left(19\right)\end{array}$

[0082] where D_nis the diffusion constant of electrons, D_pis the dispersion constant of holes, μ_nis the drift mobility of electrons, μ_pis the drift mobility of holes, G is the generation rate of electrons, and R is the recombination rate of electrons and holes.

[0083] In the impurity diffusion simulation of the first embodiment, different equations are defined for the oxide film region and the silicon region in diffusion of boron. In this device simulation, however, the definition regions of the Poisson equation and the current continuous equations are partially different.

[0084] The simulation operation of the second embodiment will be described below with reference to the flow chart shown in FIG. 2.

[0085] First, the number of equations defined at each mesh point is checked (step 204). In the second embodiment; three equations, e.g., the above-described equations (17), (18), and (19) are defined at each of the mesh points 604 on the silicon substrate side in FIG. 6. In addition, the above-described equation (17) is defined at each of the mesh points 605 on the surface of the oxide film region. Subsequently, the equations to be solved are put into a group in units of mesh points to form a coefficient matrix represented by blocks shown in FIG. 7 (step 205). In FIG. 7, a coefficient matrix having 18 rows and 18 columns in correspondence with numbers “1” to “18” assigned to the mesh points shown in FIG. 6 is obtained. In one row of the coefficient matrix shown in FIG. 7, all of the square submatrix blocks at the mesh points constituting the principal diagonals having the row numbers as the mesh point numbers shown in FIG. 6 and the square or non-square submatrix blocks at adjacent mesh points connected to the mesh point are arranged.

[0086] In the second embodiment, for the mesh point 604, the width of each block is 3, as indicated by row “4” and column “4” to row “18” and column “18”. To the contrary, for the mesh points 605 with numbers “1” to “3”, the width of each block is 1, as indicated by row “1” and column “1” to row “3” and column “3”. In each block shown in FIG. 7, a blank portion represents that a numerical value to be calculated is not stored.

[0087] For example, for row “5” in FIG. 7, which corresponds to the mesh point 604 represented by “5” in FIG. 6, a principal diagonal block of a 3×3 submatrix of the mesh point 604 with “5” is present at row “5” and column “5”. This principal diagonal block has 9 numerical values to be calculated. As shown in FIG. 6, the mesh point 604 represented by “8” is present under the mesh point 604 represented by “5”, so a 3×3 submatrix block corresponding to the mesh point 604 represented by “8” is present at row “5” and column “8” in FIG. 7. This submatrix block has five numerical values to be calculated. The mesh point 604 represented by “4” and the mesh point 604 represented by “6” are present on the left and right sides of the mesh point 604 represented by “5”, respectively, so 3×3 submatrix blocks corresponding to these mesh points are present at row “5” and column “4” and row “5” and column “6”. Each of these submatrix blocks also has five numerical values to be calculated. The mesh point 605 represented by “2” is present above the mesh point 604 represented by “5”, so a 3×1 non-square submatrix block corresponding to the mesh point 605 is present at row “5” and column “2” shown in FIG. 7.

[0088] Similarly, for row “2” corresponding to the mesh point 605 represented by “2” in FIG. 6, a principal diagonal block of a 1×1 submatrix corresponding to the mesh point 605 is present at row “2” and column “2”. Since the mesh points 605 represented by “1” and “3” are present on the left and right sides of the mesh point 605 represented by “2”, 1×1 submatrix blocks are present at row “2” and column “1” and row “2” and column “3” in FIG. 7. In addition, since the mesh point 604 represented by “5” is present under the mesh point 605 represented by “2”, a 1×3 non-square submatrix block is present at row “2” and column “5” in FIG. 7.

[0089] Finally, each submatrix is regarded as one element, and calculation for square matrices and calculation for non-square matrices are performed to perform incomplete LU-factorization in accordance with equations (1) and (2), thereby solving the equations and simulating the states at the above-described mesh points (step 206).

[0090] Conventionally, since a 3×3 submatrix is formed for all mesh points, an excess memory capacity is used, and wasteful calculation is performed. In the second embodiment, however, calculation is performed using matrices whose calculation amount is partially smaller than that of 3×3 submatrices, so wasteful calculation is omitted, and the calculation time can be shortened.

[0091] The present invention has been described with reference to the preferred embodiments, but is not limited to them.

[0092] According to the present invention, there is provided the semiconductor process device simulation method comprising a first step of dividing the surface of a semiconductor device to be processed into a mesh of predetermined size, a second step of assigning a numerical value to each mesh point of the mesh, a third step of setting equations representing a relationship among the plurality of numerical values, a fourth step of generating a coefficient matrix constituted by a principal diagonal submatrix arranged at any one of principal diagonals corresponding to each mesh point and representing a self feedback function at the mesh point, the principal diagonal submatrix having rows and columns in numbers corresponding to the number of mesh points, and a non-principal diagonal submatrix arranged on any one of a row and column passing through principal diagonal positions corresponding to the mesh point and representing an interaction between the mesh point corresponding to the principal diagonal positions and an adjacent mesh point connected to the mesh point through a mesh branch, and a fifth step of performing calculation for the submatrices while regarding each submatrix of the coefficient matrix as one element, thereby performing incomplete LU-factorization of the coefficient matrix, wherein the principal diagonal submatrix is a square matrix having rows and columns equal in number to equations set for a mesh point corresponding to the principal diagonal submatrix, the non-principal diagonal submatrix being arranged in a row direction in correspondence with each of the mesh points is a matrix having rows equal in number to equations set at a mesh point corresponding to the principal diagonal submatrix and columns equal in number to equations set at an adjacent mesh point connected to the mesh point through a mesh branch, and the non-principal diagonal submatrix being arranged in a column direction in correspondence with each of the mesh points is a matrix having columns equal in number to equations set at a mesh point corresponding to the principal diagonal submatrix and rows equal in number to equations set at an adjacent mesh point connected to the mesh point through a mesh branch.

[0093] With this arrangement, the non-principal diagonal submatrix is not always a square submatrix. When the number of equations set at an adjacent mesh point connected to a mesh point of interest through the mesh branch is smaller than that at the mesh point of interest, a non-square matrix is obtained, so the number of times of calculation becomes smaller than that when a square matrix is simply assigned. Therefore, in the present invention, less memory is required for calculation, and the simulation result can be obtained at a high speed.