Method for simulating diffusion of impurity materials in semiconductor devices

A method for simulating the deformation of regions in a semiconductor device due to oxidation. An oxidation calculation triangular mesh is deformed according to an oxidation calculation, and a diffusion calculation triangular mesh simulates the diffusion of impurity materials. Diffusion calculation control volumes are defined to each vertex of the diffusion calculation triangular mesh. The diffusion calculation triangular mesh and the diffusion calculation control volumes are deformed according to the deformation of the oxidation calculation triangular mesh. Impurity material concentrations are altered according to volume ratio of the diffusion calculation control volume after deformation to the diffusion calculation control volume before deformation. A new diffusion calculation triangular mesh and corresponding new diffusion calculation control volumes are defined for deformed shapes of regions that guarantee Delaunay partitioning. The amount of the impurity material in the diffusion calculation control volumes after deformation is fully transferred to the new diffusion calculation control volumes. A diffusion calculation is executed using the new diffusion calculation control volumes and the impurity material concentrations on the new diffusion calculation triangular mesh. The above process is repeated for a predetermined period of time. The simulation conserves the gross dose of the impurity material with a high degree of accuracy and a short calculation time.

BACKGROUND OF THE INVENTION 
The present invention relates to a method for process simulation of a 
semiconductor device, and in particular, to a process simulation method in 
which process simulation keeping conservation of the gross amount of 
impurity materials in the device is possible with short calculation time 
even in the case where the shapes of regions of the semiconductor device 
vary due to oxidation calculation etc. 
DESCRIPTION OF THE PRIOR ART 
In the field of process simulation, process simulators are used on computer 
systems for calculating and simulating manufacturing processes of 
semiconductor devices such as an oxidation process, a diffusion process, 
an ion implantation process, etc. and predicting or estimating physical 
quantities or shapes such as impurity profile of a transistor etc. For 
example, in the case where oxidation and diffusion of impurity materials 
progress in a wafer having initial shapes of parts of a semiconductor 
device right after ion implantation of the impurity materials, it is 
possible to predict time dependent variation or transition of impurity 
profile in the device and final shape of parts of the device, by 
alternately solving variation of shapes due to oxidation at each 
oxidation/diffusion time and diffusion of the impurity material (i.e. 
flux) in the atmosphere for oxidation. 
By optimizing transistors using the process simulator so that semiconductor 
devices composed of the transistors may display the best performance or 
electrical property, both cost and time needed for manufacturing the 
semiconductor devices can be considerably reduced, compared with the case 
in which a lot of LSIs etc. are actually manufactured on a trial basis. 
Furthermore, it is possible to make analyses of movement of impurity 
materials in semiconductor devices, since process simulators are capable 
of calculating internal physical quantities in semiconductor devices. In 
order to obtain internal physical quantities by process simulators, it is 
needed to solve partial differential equations such as diffusion 
continuity equation representing movement of impurity materials. However, 
it is impossible to analytically solve these kinds of partial differential 
equations. Therefore, in use of process simulators, a semiconductor device 
is divided or partitioned into a plurality of small elements, and partial 
differential equations are applied to these small elements discretely, 
As an example of a method for calculating one-dimensional impurity profile 
due to oxidation/diffusion, there is a document for reference by Morisue: 
"Designing and Manufacture Simulation of VLSI" (CMC, pages 51-62). 
And as an example of a method for analyzing two-dimensional structures, 
there is a method in which a semiconductor device is partitioned into a 
plurality of small rectangular elements and partial differential equations 
are applied to these small rectangular elements discretely, as described 
in `Process, Device, and Simulation Technique` (Sangyo Tosho, pages 
90-122) by Dan. 
Meanwhile, in order to execute analysis of a device having complicated 
shapes such as LOCOS structure or trench structure, there exists a device 
simulation method in which a semiconductor device is partitioned into a 
plurality of discrete triangular elements so that the complicated shapes 
of regions of the semiconductor device can be precisely realized in the 
simulation, as disclosed in C. S. Rafferty et al. `Iterative Methods in 
Semiconductor Device Simulation` (IEEE Trans. on ED, Vol.ED-32, No. 10, 
October 1985, pages 2018-2027). 
FIG. 1A through FIG. 1C are schematic diagrams showing an example of a 
simulation for a semiconductor device having trench structure which is 
quoted from aforementioned document of C. S. Rafferty et al. Each figure 
shows each process of a simulation for a CMOS which has trench isolated 
structure, in which FIG. 1A is showing a sectional structure of the CMOS. 
FIG. 1B is showing initial grid for the simulation, and FIG. 1C is showing 
generated working grid. By partitioning the shapes of parts of a device 
into small discrete elements using the triangular grid, the shapes of 
parts are expressed by a set of triangular elements and the shape of the 
trench structure can be precisely expressed. 
In the following, solution for the aforementioned partial differential 
equations by finite difference method using the triangular grid will be 
briefly explained. 
FIG. 2 is an enlarged schematic diagram showing part of FIG. 1B for 
explaining the finite difference method utilizing the triangular elements. 
First, concentration of an impurity material and an electric potential due 
to activated impurity material are defined on each lattice point of the 
triangular grid, i.e. on the vertexes of triangles. The impurity material 
diffuses according to a concentration gradient of the impurity material 
and an electric potential gradient (i.e. an electric field) at the point. 
The flux of the diffusing impurity material is defined on each side of the 
triangles. Here, according to theorem of Gauss, when a closed surface is 
defined, a volume integral of the impurity material concentration in the 
closed surface is equal to a surface integral of normal components of the 
impurity material flux all over the closed surface. 
The theorem of Gauss is applied to each `volume element` corresponding to 
each vertex of the triangular grid, in which each volume element has a 
polygonal shape and unit thickness in three-dimensional space. In order to 
apply the theorem of Gauss to the discrete volume element, the closed 
surface, i.e. the surface of the volume element, is needed to be defined 
perpendicular to the flux, that is , each side of the polygon has to be 
perpendicular to the flux. Therefore, the closed surface, i.e. the surface 
of the volume element, in the theorem of Gauss has to be defined so that 
the polygon is surrounded by perpendicular bisectors of the triangles, 
that is, the polygon is surrounded by sides each of which is connecting 
between circumcenters of the triangles. Here, each volume element which is 
corresponding to each vertex of the triangles is generally called `a 
control volume`. In this case, the gross amount of the impurity material 
controlled by each vertex is equal to the concentration of the impurity 
material at the vertex multiplied by the volume (the area.times.1, in two 
dimensional models) of the control volume, and the sum of the gross 
amounts corresponding to all the vertexes in the system to be analyzed is 
equal to the gross dose of the impurity material in the ion implantation. 
Incidentally, in order to realize appropriate control volumes (i.e. 
appropriate close surfaces corresponding to each vertex), there exists a 
necessary condition that the distance between circumcenters of adjacent 
triangles should not be negative, since negative cross section for 
integration occurs in the surface integral of the flux if the distance 
between circumcenters of the adjacent triangles becomes negative. If the 
necessary condition is not satisfied, physically impossible voltage spikes 
occur as shown in FIG. 3 which is quoted from aforementioned document of 
C. S. Rafferty et al. Further, if the necessary condition is not 
satisfied, overlap of more than one control volume occurs, and thus the 
aforementioned sum of the impurity material gross amounts corresponding to 
all the vertexes in the system can not keep being equal to the gross dose 
of the ion implantation. 
For satisfying the necessary condition that the distance between 
circumcenters of adjacent triangles should not be negative, it is 
necessary to partition two-dimensional regions into the triangular 
elements guaranteeing `Delaunay partitioning` which requires that no 
vertexes of other triangles should not exist in a circumscribed circle of 
a triangle. There is a known method for partitioning the region into 
triangular elements guaranteeing the Delaunay partitioning, as disclosed 
in M. S. Mock "Tetrahedral Elements and the Scharferrer-Gummel Methods" 
(Proceedings of the NASECODE IV, pages 36-47, 1985). 
Next, in the following, a conventional calculation method for simulation of 
oxidation will be described. 
As described in "Semiconductor Process.Device Simulation Technique" 
(REALIZE INC., 1990), pages 79-89, S. Isomae "Simulation of 
Two-Dimensional Oxidation", oxidation phenomenon can be considered to be 
comprised of two alternate processes: diffusion of an oxidant or an 
oxidizing agent in the oxide layer and variation of the shape of the oxide 
layer due to oxidation (for example, Si.fwdarw.SiO.sub.2). 
In calculation of time-dependent variation of shapes of regions of the 
semiconductor device, first, distribution of concentration of the oxidant 
in the oxide layer is obtained by solving the two-dimensional Laplace's 
differential equation (1) in the shape before deformation due to 
oxidization. 
EQU D.gradient..sup.2 C (x, y)=0 (1) 
Here, D: diffusion coefficient of the oxidant in the oxide layer, C(x,y): 
concentration of the oxidand at the position (x,y). 
Subsequently, variation of the shapes is calculated. Movement distance of 
each vertex of the triangular grid (triangular mesh) due to oxidation is 
obtained by solving the following equilibrium equation (2) and 
constitutive equation (3), using movement distance of oxidation interface 
according proportionally to concentration of the oxidant on the oxidation 
interface obtained by the equation (1), and displacement due to inflation 
or swelling by oxidation, as boundary conditions. 
##EQU1## 
Here, .sigma..sub.il (t) is an element of a stress tensor, X.sub.j is a 
position, .epsilon..sub.kl is an element of a deformation tensor, 
G.sub.ijkl is a relaxation function, and .delta. is a symbol denoting 
partial differentiation. 
Incidentally, when oxidation process is calculated and simulated in the 
process simulator, the shapes of the regions vary as time passes, and the 
shapes of the discrete triangular elements vary. Therefore, even if 
aforementioned Delaunay partitioning requiring that `no vertexes of other 
triangles should not exist in a circumscribed circle of a triangle` is 
satisfied at the moment before oxidation, the Delaunay partitioning begins 
to be unguaranteed as the shapes of the discrete triangular elements vary 
due to deformation by oxidation. 
Thus, when a diffusion equation for obtaining diffusion of the impurity 
material is solved subsequently to oxidation (deformation) calculation, 
another triangular partitioning guaranteeing the Delaunay partitioning has 
to be executed to the shapes of the regions after deformation. Further, 
even in the case where another triangular partitioning guaranteeing the 
Delaunay partitioning is executed to the shapes of the regions after 
deformation, control volumes alsovary according to the reconstruction of 
the triangular elements. 
Therefore, in order to keep on guaranteeing the Delaunay partitioning and 
keep on conserving the gross dose of the impurity material introduced in 
ion implantation etc., a process simulation method for properly 
distributing the impurity material as described below is proposed by the 
present inventor in Japanese Patent Laid-Open Publication No. 
9-17738(1997). 
Briefly, the method comprises the following steps in order to conserve the 
gross dose of the impurity material. 
(1) Control volumes which has been defined before oxidation are deformed 
according to variation of shapes due to oxidation obtained by oxidation 
calculation (post-deformation control volumes). 
(2) Another triangular partitioning guaranteeing the Delaunay partitioning 
is executed to the shapes after deformation. 
(3) Control volumes are redefined for the new triangular grid (mesh) (new 
control volumes). 
(4) By executing AND operation between the post-deformation control volumes 
and the new control volumes, impurity material in the post-deformation 
control volumes is fully handed over to the new control volumes with no 
loss. 
FIG. 4 is a flow chart showing the steps involved in the process simulation 
method presently proposed by the present inventor. 
First, triangular mesh (grid) for the initial shapes of regions of a 
semiconductor device to be simulated is generated by the process 
simulator, guaranteeing Delaunay partitioning. Partitioning guaranteeing 
Delaunay partitioning is executed according to a similar method to the 
method elaborated in the aforementioned document "Tetrahedral Elements and 
the Scharferrer-Gummel Method" (pre-deformation triangular mesh) (step 
S1). 
Subsequently, control volumes for the pre-deformation triangular mesh are 
defined by the process simulator. In other words, the process simulator 
defines aforementioned each `volume element` corresponding to each vertex 
of the triangular grid as a control volume for the vertex of the 
pre-deformation triangular mesh. Incidentally, the surface of the `volume 
element` (i.e. the polygon) is surrounded by sides each of which is 
connecting between circumcenters of the triangles (pre-deformation control 
volumes) (step S2). 
Subsequently, in the initial shapes of regions of the semiconductor device, 
impurity material concentrations due to introduction by ion implantation 
etc. are set by the process simulator on each vertex of the 
pre-deformation triangular mesh (step S3). 
Subsequently, unit time .DELTA.t for oxidation (deformation) and diffusion 
is added to the time counted by the process simulator (step S4). 
Subsequently, oxidation calculation is executed and deformation of the 
triangular mesh according to the oxidation calculation is executed by the 
process simulator. In other words, deformation of the pre-deformation 
triangular mesh due to oxidation is calculated and positions of each 
vertex of the pre-deformation triangular mesh are displaced 
(post-deformation triangular mesh) (step S5). 
Subsequently, the control volumes corresponding to each vertex of the 
triangular mesh are deformed by the process simulator according to the 
deformation of the triangular mesh (post-deformation control volumes) 
(step S6). 
Subsequently, volume ratios of the post-deformation control volumes to the 
pre-deformation control volumes are obtained by the process simulator and 
impurity material concentrations are set by the process simulator on each 
vertex of the post-deformation triangular mesh according to the obtained 
volume ratio and impurity material concentrations on each vertex of the 
pre-deformation triangular mesh. Here, impurity material concentration on 
a vertex on the post-deformation triangular mesh is set equal to impurity 
material concentration on the vertex on the pre-deformation triangular 
mesh multiplied by the obtained volume ratio, and thus impurity material 
concentrations on each vertex get lower according to control volumes which 
are swelling due to oxidation (step S7). 
Subsequently, after the setting of impurity material concentrations after 
deformation, another triangular mesh for deformed shapes of regions of the 
semiconductor device is newly generated by the process simulator, 
guaranteeing the Delaunay partitioning (new triangular mesh) (step S8). 
Subsequently, another control volumes corresponding to each vertex are 
defined for the new triangular mesh by the process simulator (new control 
volumes) (step S9). 
Subsequently, the amount of the impurity material in the post-deformation 
control volumes is fully transferred to the new control volumes by the 
process simulator by executing AND operation between the post-deformation 
control volumes and the new control volumes (step S10). 
Subsequently, diffusion calculation is executed by the process simulator 
using the new control volumes and the concentrations of the impurity 
material on vertexes of the new triangular mesh, i.e. diffusion equation 
of the impurity material is solved (step S11). 
Subsequently, it is judged whether or not predetermined time has been 
counted. If NO, process is returned to the step S4 to repeat the step S4 
through the step S11 until the predetermined time is counted by the 
process simulator. If YES, the process simulation is ended (step S12). 
FIG. 5 is a flow chart showing the process of the step S10 of FIG. 4 in 
which the amount of the impurity material in the post-deformation control 
volumes is fully transferred to the new control volumes by executing AND 
operation. In the following, the steps involved in the process of the step 
S10 will be described. 
First, a new control volume ICVnew is arbitrarily chosen (step S21). 
Subsequently, a post-deformation control volume ICVold corresponding to 
the chosen new control volume ICVnew is chosen (step S22). 
Subsequently, AND operation is executed between the ICVnew and the ICVold, 
and the area S.sub.AND of the overlap between ICVnew and the ICVold is 
obtained (step S23). Subsequently, the amount of the impurity material 
D.sub.AND in the overlap is obtained by multiplying the S.sub.AND and the 
impurity material concentration set on the post-deformation triangular 
mesh in the step S7 of FIG. 4 (step S24). Subsequently, obtained D.sub.AND 
is added to the amount of impurity material Dnew in the new control volume 
ICVnew (step S25). 
Subsequently, it is judged whether or not the step S23 through the step S25 
has been executed to all of post-deformation control volumes overlapping 
with the new control volume (step S26). If NO, next post-deformation 
control volumes overlapping with the new control volume is chosen, and 
process is returned to the step S23 in order to repeat the step S23 
through the step S25 (step S29). If YES, impurity material concentration 
on a vertex corresponding to the new control volume ICVnew is obtained by 
dividing the resultant Dnew by the area of the new control volume ICVnew 
(step S27). 
Subsequently, it is judged whether or not the aforementioned process has 
been executed to all of new control volumes (step S28). If NO, next new 
control volume ICVnew is chosen, and process is returned to the step S22 
in order to repeat the step S22 through the step S29 (step S30). 
As described above, according to the process simulation method presently 
proposed by the present inventor, the pre-deformation control volumes 
corresponding to the pre-deformation triangular mesh is deformed to the 
post-deformation control volumes corresponding to the post-deformation 
triangular mesh according to deformation of the triangular mesh obtained 
in the oxidation calculation, impurity material concentrations on vertexes 
of the post-deformation triangular mesh are properly set according to 
impurity material concentrations on vertexes of the pre-deformation 
triangular mesh and the obtained volume ratio due to deformation, the new 
triangular mesh is generated guaranteeing Delaunay partitioning and the 
new control volumes are defined for deformed shapes of regions of the 
semiconductor device, and the diffusion calculation is executed after the 
amount of the impurity material in the post-deformation control volumes is 
fully and properly transferred to the new control volumes by executing AND 
operation between post-deformation control volumes and new control 
volumes. 
Therefore, by the process simulation method presently proposed by the 
present inventor, process simulation keeping on guaranteeing the Delaunay 
partitioning and keeping conservation of the gross dose of the impurity 
material due to introduction such as ion implantation is realized, even in 
the case where shapes of regions of the semiconductor device vary as time 
passes due to oxidation. 
However, the aforementioned method proposed by the present inventor has a 
drawback to need considerable calculation time, since the same triangular 
mesh is used for both oxidation calculation and diffusion calculation, and 
thus fine mesh is needed both in the silicon oxide layer and in the 
silicon layer in order to execute calculation with high precision. 
For both high precision calculation and short calculation time, the 
triangular mesh for oxidation calculation is needed to be fine in and 
around the silicon oxide layer and coarse in the silicon layer. On the 
other hand, the triangular mesh for diffusion calculation is needed to be 
finer in the silicon layer and coarser in the silicon oxide layer compared 
to the triangular mesh for oxidation calculation. However, in the 
aforementioned proposed method, triangular mesh after diffusion 
calculation is utilized as triangular mesh for the next oxidation 
calculation, and alteration of density of the triangular mesh is 
impossible. 
It is of course possible to generate triangular mesh for oxidation 
calculation after diffusion calculation. However, it is necessary to 
transfer and properly distribute the impurity material from the triangular 
mesh for diffusion calculation into the triangular mesh for oxidation 
calculation. In this case, when the impurity material is distributed from 
the triangular mesh for diffusion calculation in which the silicon layer 
has fine mesh into the triangular mesh for oxidation calculation in which 
the silicon layer has coarse mesh, detailed information of impurity 
profile in the silicon layer is lost. 
SUMMARY OF THE INVENTION 
It is therefore the primary object of the present invention to provide a 
process simulation method by which process simulation keeping on 
guaranteeing the Delaunay partitioning and keeping on conserving the gross 
dose of the impurity material due to introduction such as ion implantation 
is realized even in the case where shapes of regions of the semiconductor 
device vary as time passes due to oxidation, with both high precision 
calculation and short calculation time. 
In accordance with the present invention, there is provided a process 
simulation method comprising the following eleven steps. In the first 
step, pre-deformation diffusion calculation triangular mesh is generated 
for initial shapes of regions of a semiconductor device guaranteeing 
Delaunay partitioning. In the second step, pre-deformation diffusion 
calculation control volumes are defined for the pre-deformation diffusion 
calculation triangular mesh. In the third step, impurity material 
concentrations due to introduction by ion implantation etc. are set on 
each vertex of the pre-deformation diffusion calculation triangular mesh 
in the initial shapes of regions of the semiconductor device. 
In the fourth step, oxidation calculation triangular mesh is generated for 
the shapes of regions of the semiconductor device before oxidation. In the 
fifth step, oxidation calculation is executed and the oxidation 
calculation triangular mesh is deformed according to the oxidation 
calculation. In the sixth step, pre-deformation diffusion calculation 
triangular mesh and pre-deformation diffusion calculation control volumes 
are deformed to post-deformation diffusion calculation triangular mesh and 
post-deformation diffusion calculation control volumes respectively 
according to the deformation of the oxidation calculation triangular mesh 
in the fifth step. In the seventh step, impurity material concentrations 
on each vertex of the post-deformation diffusion calculation triangular 
mesh are set according to impurity material concentrations on each vertex 
of the pre-deformation diffusion calculation triangular mesh and volume 
ratio of the post-deformation diffusion calculation control volume to the 
pre-deformation diffusion calculation control volume. In the eighth step, 
new diffusion calculation triangular mesh is generated for deformed shapes 
of regions of the semiconductor device guaranteeing Delaunay partitioning. 
In the ninth step, new diffusion calculation control volumes are defined 
for the new diffusion calculation triangular mesh. In the tenth step, the 
amount of the impurity material in the post-deformation diffusion 
calculation control volumes is fully transferred to the new diffusion 
calculation control volumes by executing AND operation between the 
post-deformation diffusion calculation control volumes and the new 
diffusion calculation control volumes. And in the eleventh step, diffusion 
calculation is executed using the new diffusion calculation control 
volumes and the concentrations of the impurity material on vertexes of the 
new diffusion calculation triangular mesh. The process between the fourth 
step and the eleventh step for unit time.DELTA.t is repeated until 
predetermined time is counted by the process simulator. 
Preferably, the sixth step mentioned above includes the following eight 
steps. In the first step, a position P.sub.diff of a vertex or a 
circumcenter of a triangle composing the pre-deformation diffusion 
calculation triangular mesh is arbitrarily chosen. In the second step, a 
triangle T.sub.ox composing the oxidation calculation triangular mesh 
before deformation is chosen. In the third step, it is judged whether or 
not the chosen triangle L.sub.ox includes the chosen position P.sub.diff. 
In the fourth step, if the judgment in the third step is NO, next triangle 
T.sub.ox composing the oxidation calculation triangular mesh before 
deformation is chosen and process is returned to the third step. In the 
fifth step, if the judgment in the third step is YES, the position 
P.sub.diff is moved using movement of three vertexes of the triangle 
T.sub.ox of the oxidation calculation triangular mesh by deformation due 
to oxidation. In the sixth step, it is judged whether or not all of the 
vertexes and the circumcenters of the pre-deformation diffusion 
calculation triangular mesh have been moved. In the seventh step, if the 
judgment in the sixth step is NO, next position P.sub.diff of a vertex or 
a circumcenter of a triangle composing the pre-deformation diffusion 
calculation triangular mesh is chosen and process is returned to the 
second step. And in the eight step, if the judgment in the sixth step is 
YES, post-deformation diffusion calculation control volumes are defined 
for the resultant post-deformation diffusion calculation triangular mesh. 
Preferably, in the fifth step mentioned above, in the case where positions 
of the three vertexes of the triangle T.sub.ox in the oxidation 
calculation triangular mesh: 
##EQU2## 
by deformation due to oxidation, the position P.sub.diff of the vertex or 
the circumcenter in the diffusion calculation triangular mesh: P.sub.diff 
(X.sub.diff, Y.sub.diff) is moved into P.sub.diff.sup.move 
(X.sub.diff.sup.move, Y.sub.diff.sup.move) as the following equations. 
##EQU3## 
Preferably, the aforementioned tenth step includes the following eleven 
steps. In the first step, a new diffusion calculation control volume is 
chosen. In the second step, a post-deformation diffusion calculation 
control volume corresponding to the chosen new diffusion calculation 
control volume is chosen. In the third step, AND operation is executed 
between the chosen new diffusion calculation control volume and the chosen 
post-deformation diffusion calculation control volume, and the area S 
.sub.AND of the overlap between the two control volumes is obtained. In 
the fourth step, the amount D.sub.AND of the impurity material in the 
overlap is obtained by multiplying the S.sub.AND and the impurity material 
concentration set on the post-deformation diffusion calculation triangular 
mesh. In the fifth step, obtained D.sub.AND is added to the amount of the 
impurity material in the new diffusion calculation control volume. In the 
sixth step, it is judged whether or not the third step through the fifth 
step have been executed to all post-deformation diffusion calculation 
control volumes overlapping with the new diffusion calculation control 
volume. In the seventh step, if the judgment in the sixth step is NO, next 
post-deformation diffusion calculation control volume overlapping with the 
new diffusion calculation control volume is chosen and process is returned 
to the third step. In the eighth step, if the judgment in the sixth step 
is YES, impurity material concentration on a vertex corresponding to the 
new diffusion calculation control volume is obtained by dividing the 
resultant amount of the impurity material in the new diffusion calculation 
control volume by the area of the new control volume. In the ninth step, 
it is judged whether or not the second step through the eighth step have 
been executed to all new diffusion calculation control volumes. In the 
tenth step, if the judgment in the ninth step is NO, next new diffusion 
calculation control volume is chosen and process is returned to the second 
step. And in the eleventh step, if the judgment in the ninth step is YES, 
the process for the aforementioned tenth step is finished. 
In accordance with another aspect of the present invention, there is 
provided a process simulation method comprising the following twelve 
steps. In the first step, pre-deformation diffusion calculation triangular 
mesh is generated for initial shapes of regions of a semiconductor device 
guaranteeing Delaunay partitioning. In the second step, pre-deformation 
diffusion calculation control volumes are defined for the pre-deformation 
diffusion calculation triangular mesh. In the third step, impurity 
material concentrations due to introduction by ion implantation etc. are 
set on each vertex of the pre-deformation diffusion calculation triangular 
mesh in the initial shapes of regions of the semiconductor device. In the 
fourth step, vertexes of the oxidation calculation triangular mesh are 
added to points where stress gradient is large. In the fifth step, 
oxidation calculation triangular mesh for the shapes of regions of the 
semiconductor device before oxidation is generated using the vertexes 
which has been added in the fourth step. In the sixth step, oxidation 
calculation is executed and the oxidation calculation triangular mesh is 
deformed according to the oxidation calculation. In the seventh step, 
pre-deformation diffusion calculation triangular mesh and pre-deformation 
diffusion calculation control volumes are deformed to post-deformation 
diffusion calculation triangular mesh and post-deformation diffusion 
calculation control volumes respectively according to the deformation of 
the oxidation calculation triangular mesh in the sixth step. In the eighth 
step, impurity material concentrations on each vertex of the 
post-deformation diffusion calculation triangular mesh are set according 
to impurity material concentrations on each vertex of the pre-deformation 
diffusion calculation triangular mesh and volume ratio of the 
post-deformation diffusion calculation control volume to the 
pre-deformation diffusion calculation control volume. In the ninth step, 
new diffusion calculation triangular mesh is generated for deformed shapes 
of regions of the semiconductor device guaranteeing Delaunay partitioning. 
In the tenth step, new diffusion calculation control volumes are defined 
for the new diffusion calculation triangular mesh. In the eleventh step, 
the amount of the impurity material in the post-deformation diffusion 
calculation control volumes is fully transferred to the new diffusion 
calculation control volumes by executing AND operation between the 
post-deformation diffusion calculation control volumes and the new 
diffusion calculation control volumes. And in the twelfth step, diffusion 
calculation is executed using the new diffusion calculation control 
volumes and the concentrations of the impurity material on vertexes of the 
new diffusion calculation triangular mesh. The process between the fourth 
step and the twelfth step for unit time.DELTA.t is repeated until 
predetermined time is counted by the process simulator. 
Preferably, each of the vertexes added in the step fourth step is added to 
a midpoint of adjacent two vertexes of the oxidation calculation 
triangular mesh if the stress ratio between the adjacent two vertexes is 
larger than a predetermined value. 
Preferably, the predetermined value is two. 
In accordance with another aspect of the present invention, there is 
provided a process simulation method comprising the following twelve 
steps. In the first step, pre-deformation diffusion calculation triangular 
mesh is generated for initial shapes of regions of a semiconductor device 
guaranteeing Delaunay partitioning. In the second step, pre-deformation 
diffusion calculation control volumes are defined for the pre-deformation 
diffusion calculation triangular mesh. In the third step, impurity 
material concentrations due to introduction by ion implantation etc. are 
set on each vertex of the pre-deformation diffusion calculation triangular 
mesh in the initial shapes of regions of the semiconductor device. In the 
fourth step, oxidation calculation triangular mesh is generated for the 
shapes of regions of the semiconductor device before oxidation. In the 
fifth step, oxidation calculation is executed and the oxidation 
calculation triangular mesh is deformed according to the oxidation 
calculation. In the sixth step, pre-deformation diffusion calculation 
triangular mesh and pre-deformation diffusion calculation control volumes 
are deformed to post-deformation diffusion calculation triangular mesh and 
post-deformation diffusion calculation control volumes respectively 
according to the deformation of the oxidation calculation triangular mesh 
in the fifth step. In the seventh step, impurity material concentrations 
on each vertex of the post-deformation diffusion calculation triangular 
mesh are set according to impurity material concentrations on each vertex 
of the pre-deformation diffusion calculation triangular mesh and volume 
ratio of the post-deformation diffusion calculation control volume to the 
pre-deformation diffusion calculation control volume. In the eighth step, 
vertexes of the diffusion calculation triangular mesh are added to points 
where concentration gradient of the impurity material is large. In the 
ninth step, new diffusion calculation triangular mesh for deformed shapes 
of regions of the semiconductor device is generated using the vertexes 
which has been added in the eighth step guaranteeing Delaunay 
partitioning. In the tenth step, new diffusion calculation control volumes 
are defined for the new diffusion calculation triangular mesh. In the 
eleventh step, the amount of the impurity material in the post-deformation 
diffusion calculation control volumes is fully transferred to the new 
diffusion calculation control volumes by executing AND operation between 
the post-deformation diffusion calculation control volumes and the new 
diffusion calculation control volumes. And in the twelfth step, diffusion 
calculation is executed using the new diffusion calculation control 
volumes and the concentrations of the impurity material on vertexes of the 
new diffusion calculation triangular mesh. The process between the fourth 
step and the twelfth step for unit time.DELTA.t is repeated until 
predetermined time is counted by the process simulator. 
Preferably, each of the vertexes added in the eighth step is added to a 
midpoint of adjacent two vertexes of the diffusion calculation triangular 
mesh if the concentration ratio of the impurity material between the 
adjacent two vertexes is larger than a predetermined value. 
Preferably, the predetermined value is ten. 
In accordance with another aspect of the present invention, there is 
provided a process simulation method using oxidation calculation 
triangular mesh for simulating deformation of regions of a semiconductor 
device due to oxidation, diffusion calculation triangular mesh for 
simulating diffusion of impurity materials introduced by ion implantation 
etc. in the semiconductor device, and diffusion calculation control 
volumes defined to each vertex of the diffusion calculation triangular 
mesh, comprising the following eight steps. In the first step, the 
diffusion calculation triangular mesh is generated for shapes of regions 
of the semiconductor device before deformation due to oxidation 
guaranteeing Delaunay partitioning, the oxidation calculation triangular 
mesh is generated for the shapes before deformation due to oxidation, and 
the diffusion calculation control volumes are defined to each vertex of 
the diffusion calculation triangular mesh. In the second step, impurity 
material concentrations are set on each vertex of the diffusion 
calculation triangular mesh in the initial shapes of regions of the 
semiconductor device. In the third step, oxidation calculation is executed 
and the oxidation calculation triangular mesh is deformed according to the 
oxidation calculation. In the fourth step, the diffusion calculation 
triangular mesh and the diffusion calculation control volumes are deformed 
according to the deformation of the oxidation calculation triangular mesh 
in the third step. In the fifth step, the impurity material concentrations 
on each vertex of the diffusion calculation triangular mesh are altered 
according to volume ratio of the diffusion calculation control volume 
after deformation to the diffusion calculation control volume before 
deformation. In the sixth step, new diffusion calculation triangular mesh 
and corresponding new diffusion calculation control volumes are defined 
for deformed shapes of regions of the semiconductor device guaranteeing 
Delaunay partitioning. In the seventh step, the amount of the impurity 
material in the diffusion calculation control volumes after deformation is 
fully transferred to the new diffusion calculation control volumes by 
executing AND operation between the diffusion calculation control volumes 
after deformation and the new diffusion calculation control volumes. And 
in the eighth step, diffusion calculation is executed using the new 
diffusion calculation control volumes and the concentrations of the 
impurity material on vertexes of the new diffusion calculation triangular 
mesh. The process between the third step and the eighth step for unit time 
.DELTA.t is repeated until predetermined time is counted by the process 
simulator.

DESCRIPTION OF THE PREFERRED EMBODIMENTS 
Referring now to the drawings, a description will be given in detail of 
preferred embodiments in accordance with the present invention. 
FIG. 6 is a flow chart showing the steps in the process simulation method 
according to an embodiment of the present invention. FIG. 7A through FIG. 
7E are schematic diagrams showing variation of triangular mesh according 
to the embodiment. The process simulation method shown in FIG. 6 is 
executed by a process simulator realized by a computer system including a 
CPU, a storage means (ROM, RAM, HDD), a software for process simulation, 
etc. 
First, triangular mesh for diffusion calculation is generated by the 
process simulator for initial shapes of regions of a semiconductor device. 
guaranteeing Delaunay partitioning. Partitioning guaranteeing Delaunay 
partitioning is executed according to a similar method to the method 
elaborated in the aforementioned document "Tetrahedral Elements and the 
Scharferrer-Gummel Method" and is not elaborated here (pre-deformation 
diffusion calculation triangular mesh) (step S31). 
FIG. 7A shows the triangular mesh for diffusion calculation generated for 
the initial shapes before oxidation, in which (1a) denotes a silicon 
nitride layer, (1b) denotes a silicon oxide layer, and (1c) denotes a 
silicon layer. In the mesh for diffusion calculation, relatively fine mesh 
is defined in the silicon layer. 
Subsequently, control volumes for the pre-deformation diffusion calculation 
triangular mesh are defined by the process simulator. In other words, the 
process simulator defines aforementioned each `volume element` 
corresponding to each vertex of the triangular grid as a control volume 
for the vertex of the pre-deformation diffusion calculation triangular 
mesh. Incidentally, the surface of the `volume element` (i.e. the polygon) 
is surrounded by sides each of which is connecting between circumcenters 
of the triangles (pre-deformation diffusion calculation control volumes) 
(step S32). 
Subsequently, in the initial shapes of regions of the semiconductor device, 
impurity material concentrations due to introduction by ion implantation 
etc. are set by the process simulator on each vertex of the 
pre-deformation diffusion calculation triangular mesh (step S33). 
Subsequently, unit time .DELTA.t for oxidation (deformation) and diffusion 
is added to the time counted by the process simulator (step S34). 
Subsequently, triangular mesh for oxidation calculation is generated by the 
process simulator for the shapes of regions before oxidation, guaranteeing 
Delaunay partitioning (pre-deformation oxidation calculation triangular 
mesh). Incidentally, control volumes are defined for diffusion calculation 
only, and definition of control volumes corresponding to the 
pre-deformation oxidation calculation triangular mesh is not executed by 
the process simulator (step S35). 
FIG. 7B shows the triangular mesh for oxidation calculation generated for 
the initial shapes before oxidation. 
In the mesh for oxidation calculation, fine mesh is defined in the silicon 
oxide layer and the silicon nitride layer in which deformation due to 
oxidation occurs. 
Subsequently, oxidation calculation is executed and deformation of the 
oxidation calculation triangular mesh according to the oxidation 
calculation is executed by the process simulator. In other words, 
deformation of the pre-deformation oxidation calculation triangular mesh 
due to oxidation is calculated and positions of each vertex of the 
pre-deformation oxidation calculation triangular mesh are displaced 
(post-deformation oxidation calculation triangular mesh) (step S36). 
Oxidation calculation in the step S36 can be executed as described in the 
description of the prior art according to the method elaborated in the 
aforementioned document "Semiconductor Process. Device Simulation 
Technique" (REALIZE INC., 1990), pages 79-89, S. Isomae "Simulation of 
Two-Dimensional Oxidation". As mentioned above, oxidation phenomenon can 
be treated as two alternate processes: diffusion of an oxidant in the 
oxide layer and variation of the shape of the oxide layer due to 
oxidation. By solving Laplace's differential equation (1) describing 
diffusion of the oxidant in oxide (SiO.sub.2) layer during unit time 
.DELTA.t, concentration of the oxidant on the oxidation interface 
(SiO.sub.2 /Si interface) is obtained. 
EQU D.gradient..sup.2 C(x, y)=0 (1) 
Here, D: diffusion coefficient of the oxidant in the oxide layer, C(x,y): 
concentration of the oxidant at the position (x,y). 
Deformation of the oxide layer is obtained by solving the equilibrium 
equation (2) and the constitutive equation (3), using movement distance of 
oxidation interface according proportionally to the concentration of the 
oxidant on the oxidation interface obtained by the equation (1), and 
displacement due to inflation or swelling by oxidation, as boundary 
conditions. 
##EQU4## 
Here, .sigma..sub.ij (t): an element of a stress tensor, x.sub.j : a 
position, .epsilon..sub.kl : an element of a deformation tensor, 
G.sub.ijkl : a relaxation function, and .delta.: a symbol denoting partial 
differentiation. 
FIG. 7C shows the post-deformation oxidation calculation triangular mesh 
and deformed shapes of regions of the semiconductor device. 
Subsequently, the pre-deformation diffusion calculation control volumes 
corresponding to each vertex of the pre-deformation diffusion calculation 
triangular mesh are deformed according to the deformation of the oxidation 
calculation triangular mesh (post-deformation diffusion calculation 
control volumes). In other words, each vertexes on the pre-deformation 
diffusion calculation triangular mesh are moved according to the 
deformation of the oxidation calculation triangular mesh to form 
post-deformation diffusion calculation triangular mesh, and then the 
post-deformation diffusion calculation control volumes are defined for 
each vertexes of the post-deformation diffusion calculation triangular 
mesh (step S37). 
FIG. 7D shows the post-deformation diffusion calculation triangular mesh in 
which positions of vertexes are moved according to the deformation of the 
oxidation calculation triangular mesh. Here, the Delaunay partitioning is 
not guaranteed in the post-deformation diffusion calculation triangular 
mesh since the mesh has been deformed. 
Subsequently, volume ratios of the post-deformation diffusion calculation 
control volumes to the pre-deformation diffusion calculation control 
volumes are obtained by the process simulator and impurity material 
concentrations are set on each vertex of the post-deformation diffusion 
calculation triangular mesh according to the obtained volume ratio and 
impurity material concentrations on each vertex of the pre-deformation 
diffusion calculation triangular mesh. Here, impurity material 
concentration on a vertex on the post-deformation diffusion calculation 
triangular mesh is set equal to impurity material concentration on the 
vertex on the pre-deformation diffusion calculation triangular mesh 
multiplied by the obtained volume ratio (step S38). 
Subsequently, after the setting of impurity material concentrations after 
deformation, another diffusion calculation triangular mesh for deformed 
shapes of regions of the semiconductor device is newly generated by the 
process simulator, guaranteeing the Delaunay partitioning (new diffusion 
calculation triangular mesh) (step S39). 
Incidentally, in order to generate the new diffusion calculation triangular 
mesh in the step 39, it is possible to use the same vertexes as those of 
the post-deformation diffusion calculation triangular mesh and alter the 
triangular elements only, to guarantee the Delaunay partitioning requiring 
that no vertexes of other triangles should not exist in a circumscribed 
circle of a triangle. However, of course it is as well possible to use new 
vertexes to generate the new diffusion calculation triangular mesh. 
FIG. 7E shows the new diffusion calculation triangular mesh generated for 
the deformed shapes of regions of the semiconductor device, in which fine 
mesh is defined in the silicon layer. 
Subsequently, another control volumes corresponding to each vertex are 
defined for the new diffusion calculation triangular mesh by the process 
simulator (new diffusion calculation control volumes) (step S40). 
Subsequently, the amount of the impurity material in the post-deformation 
diffusion calculation control volumes is fully transferred to the new 
diffusion calculation control volumes by the process simulator by 
executing AND operation between the post-deformation diffusion calculation 
control volumes and the new diffusion calculation control volumes (step 
S41). 
Subsequently, diffusion calculation is executed by the process simulator 
using the new diffusion calculation control volumes and the concentrations 
of the impurity material on vertexes of the new diffusion calculation 
triangular mesh, i.e. diffusion equation of the impurity material is 
solved (step S42). 
Subsequently, it is judged whether or not predetermined time has been 
counted. If NO, process is returned to the step S34 to repeat the step S34 
through the step S42 until the predetermined time is counted by the 
process simulator. Incidentally, in the next loop after the process is 
returned to the step S34, the new diffusion calculation triangular mesh in 
the previous loop is used as pre-deformation diffusion calculation 
triangular mesh in the step S37 of the next loop. If YES, the process 
simulation is ended (step S43). 
FIG. 8 is a flow chart showing the process of defining the post-deformation 
diffusion calculation control volumes in the step S37 of FIG. 6. In the 
following, the steps involved in the process of the step S37 will be 
described. 
First, a position P.sub.diff of a vertex or a circumcenter of a triangle 
composing the pre-deformation diffusion calculation triangular mesh is 
arbitrarily chosen (step S51). 
Subsequently, a triangle T .sub.ox composing the pre-deformation oxidation 
calculation triangular mesh is arbitrarily chosen (step S52). 
Subsequently, it is judged whether or not the chosen triangle L.sub.ox 
includes the chosen position P.sub.diff (step S53). If NO, next triangle T 
.sub.ox composing the pre-deformation oxidation calculation triangular 
mesh is arbitrarily chosen and process is returned to the step S53 (step 
S54). If YES, the position P.sub.diff is moved using movement of three 
vertexes of the triangle T.sub.ox of the pre-deformation oxidation 
calculation triangular mesh (P.sub.ox 1, P.sub.ox 2, P.sub.ox 3) (step 
S55). 
Concretely, in the case where positions of three vertexes of the triangle 
T.sub.ox in the pre-deformation oxidation calculation triangular mesh: P 
.sub.ox 1(X1.sub.ox,Y1.sub.ox), P.sub.ox 2(X2.sub.ox, Y2.sub.ox), P.sub.ox 
3(X3.sub.ox, Y3.sub.ox) moved into (X1.sub.ox.sup.move, 
Y1.sub.ox.sup.mvoe), (X2.sub.ox.sup.move, Y2.sub.ox.sup.move), 
(X3.sub.ox.sup.move, Y3.sub.ox.sup.move) by deformation due to oxidation, 
the position P.sub.diff of the vertex or the circumcenter in the diffusion 
calculation triangular mesh included in the T.sub.ox : P.sub.diff 
(X.sub.diff, Y.sub.diff) moves into P.sub.diff.sup.move 
(X.sub.diff.sup.move, Y.sub.diff.sup.move) as the following equations (4) 
and (5). 
##EQU5## 
Here, the S and T in the equations (4) and (5) is expressed by the 
following equations (6) and (7). 
##EQU6## 
Subsequently, it is judged whether or not all of the vertexes and the 
circumcenters of the pre-deformation diffusion calculation triangular mesh 
have been moved (step S56). If NO, next position P.sub.diff of a vertex or 
a circumcenter of a triangle composing the pre-deformation diffusion 
calculation triangular mesh is arbitrarily chosen and process is returned 
to the step S52 (step S57). If YES, deformation of the pre-deformation 
diffusion calculation triangular mesh into the post-deformation diffusion 
calculation triangular mesh according to the deformation of the oxidation 
calculation triangular mesh is completed, and then post-deformation 
diffusion calculation control volumes are defined according to the 
post-deformation diffusion calculation triangular mesh, and the step S37 
in FIG. 6 is completed (step S58). 
Incidentally, in the process of the step S37 of FIG. 6 described in FIG. 8, 
all of the triangular elements T.sub.oxS composing the pre-deformation 
oxidation calculation triangular mesh were referred to in order to judge 
which one of the triangular elements T.sub.oxS includes the vertexes and 
circumcenters of the pre-deformation diffusion calculation triangular 
mesh, and considerably long calculation time is needed since the scale of 
necessary algorithm is as large as the order of 0(n.sup.2) (n: the number 
of vertexes). However, in practical calculation, it is possible to execute 
high speed calculation using some methods such as well-known hash table 
method or hashing method. When the hash table method is applied to the 
step S37, the regions of the semiconductor device is previously 
partitioned into coarse rectangular elements and each triangular element 
T.sub.ox is previously registered corresponding to one of the coarse 
rectangular elements, and the coarse rectangular elements are used for 
reference, and thus calculation time is considerably shortened. 
FIG. 9 is a flow chart showing the process of the step S41 of FIG. 6 in 
which the amount of the impurity material in the post-deformation 
diffusion calculation control volumes is fully transferred to the new 
diffusion calculation control volumes by executing AND operation. In the 
following, the steps involved in the process of the step S41 will be 
described. 
First, a new diffusion calculation control volume lCVnew is arbitrarily 
chosen (step S61). Subsequently, a post-deformation diffusion calculation 
control volume ICVold corresponding to the chosen new diffusion 
calculation control volume ICVnew is chosen (step S62). Subsequently, AND 
operation is executed between the ICVnew and the ICVold, and the area 
S.sub.AND of the overlap between ICVnew and the ICVold is obtained (step 
S63). 
Subsequently, the amount of the impurity material D.sub.AND in the overlap 
is obtained by multiplying the S.sub.AND and the impurity material 
concentration set on the post-deformation diffusion calculation triangular 
mesh in the step S38 of FIG. 6 (step S64). Subsequently, obtained 
D.sub.AND is added to the amount of impurity material Dnew in the new 
diffusion calculation control volume ICVnew (step S65). 
Subsequently, it is judged whether or not the step S63 through the step S65 
have been executed to all post-deformation diffusion calculation control 
volumes overlapping with the new diffusion calculation control volume 
(step S66). If NO, next post-deformation diffusion calculation control 
volume overlapping with the new diffusion calculation control volume is 
chosen, and process is returned to the step S63 in order to repeat the 
step S63 through the step S65 (step S69). 
If YES, impurity material concentration on a vertex corresponding to the 
new diffusion calculation control volume ICVnew is obtained by dividing 
the resultant Dnew by the area of the new control volume ICVnew (step 
S67). 
Subsequently, it is judged whether or not the aforementioned process has 
been executed to all new diffusion calculation control volumes (step S68). 
If NO, next new diffusion calculation control volume ICVnew is chosen, and 
process is returned to the step S62 in order to repeat the step S62 
through the step S69 (step S70). If YES, the process of the step S41 in 
FIG. 6 is completed. 
FIG. 10 is a schematic diagram for explaining variation of the diffusion 
calculation triangular mesh and the diffusion calculation control volumes 
according to the present invention, in which A is pre-deformation 
diffusion calculation triangular mesh generated for the initial shape of a 
region of a semiconductor device guaranteeing Delaunay partitioning (step 
S31 of FIG. 6), B is a pre-deformation diffusion calculation control 
volume defined to a vertex of the pre-deformation diffusion calculation 
triangular mesh (step S32 of FIG. 6), C is post-deformation diffusion 
calculation triangular mesh which has been deformed according to 
deformation of the oxidation calculation triangular mesh due to oxidation 
(step S37 of FIG. 6), D is a post-deformation diffusion calculation 
control volume ICVold defined to the vertex of the post-deformation 
diffusion calculation triangular mesh (step S37 of FIG. 6), E is new 
diffusion calculation triangular mesh generated for deformed shape of the 
region of the semiconductor device guaranteeing Delaunay partitioning 
(step S39 of FIG. 6). Incidentally, although the same vertexes as those of 
the post-deformation diffusion calculation triangular mesh C are used in 
the new diffusion calculation triangular mesh E of FIG. 10 and only 
triangular elements are altered to guarantee the Delaunay partitioning 
requiring that no vertexes of other triangles should not exist in a 
circumscribed circle of a triangle, it is as well possible to generate the 
new diffusion calculation triangular mesh E using new vertexes. F is a new 
diffusion calculation control volume ICVnew defined to a vertex of the new 
diffusion calculation triangular mesh (step S40 of FIG. 6), and G is the 
overlap between the new diffusion calculation control volume ICVnew and 
the post-deformation diffusion calculation control volume ICVold. The AND 
operation in the step S63 of FIG. 9 is executed between the ICVnew and the 
ICVold, and the area S.sub.AND of the overlap is obtained (step S41 of 
FIG. 6, step S63 of FIG. 9). 
As described above, according to the embodiment, as the process simulation 
method presently proposed by the present inventor, process simulation 
keeping on guaranteeing the Delaunay partitioning and keeping conservation 
of the gross dose of the impurity material due to introduction such as ion 
implantation is realized, even in the case where shapes of regions of the 
semiconductor device vary as time passes due to oxidation. 
Furthermore, according to this embodiment, use of two triangular mesh, i.e. 
triangular mesh suitable for oxidation calculation and triangular mesh 
suitable for diffusion calculation, is made possible. Fine mesh can be 
defined in and around the silicon layer in the triangular mesh for 
oxidation calculation, and coarse mesh can be defined in and around the 
silicon layer in the triangular mesh for oxidation calculation, and thus 
the number of vertexes n in each triangular mesh is reduced to 
approximately 50%. The scale of necessary algorithm for oxidation 
calculation and diffusion calculation is expressed as 0(n). Therefore, 
calculation time for the process simulation can be reduced to 
approximately 50% in comparison with the case where single triangular mesh 
is used. 
FIG. 11 is a flow chart showing the steps in the process simulation method 
according to another embodiment of the present invention. 
First, pre-deformation diffusion calculation triangular mesh is generated 
by the process simulator for initial shapes of regions of a semiconductor 
device, guaranteeing Delaunay partitioning (step S71). 
Subsequently, pre-deformation diffusion calculation control volumes for the 
pre-deformation diffusion calculation triangular mesh are defined by the 
process simulator (step S72). 
Subsequently, impurity material concentrations due to introduction by ion 
implantation etc. are set by the process simulator on each vertex of the 
pre-deformation diffusion calculation triangular mesh in the initial 
shapes of regions of the semiconductor device (step S73). 
Subsequently, unit time .DELTA.t for oxidation (deformation) and diffusion 
is added to the time counted by the process simulator (step S74). 
Subsequently, vertexes of the oxidation calculation triangular mesh are 
added to points where stress gradient is large (step S75). 
Subsequently, pre-deformation oxidation calculation triangular mesh is 
generated by the process simulator for the shapes of regions before 
oxidation using the vertexes which has been added in the step S75, 
guaranteeing Delaunay partitioning (step S76). 
Subsequently, oxidation calculation is executed and deformation of the 
oxidation calculation triangular mesh according to the oxidation 
calculation is executed by the process simulator (post-deformation 
oxidation calculation triangular mesh) (step S77). 
Subsequently, the pre-deformation diffusion calculation control volumes 
corresponding to each vertex of the pre-deformation diffusion calculation 
triangular mesh are deformed according to the deformation of the oxidation 
calculation triangular mesh (post-deformation diffusion calculation 
control volumes) (step S78). 
Subsequently, volume ratios of the post-deformation diffusion calculation 
control volumes to the pre-deformation diffusion calculation control 
volumes are obtained by the process simulator and impurity material 
concentrations are set on each vertex of the post-deformation diffusion 
calculation triangular mesh according to the obtained volume ratio and 
impurity material concentrations on each vertex of the pre-deformation 
diffusion calculation triangular mesh (step S79). 
Subsequently, vertexes of the diffusion calculation triangular mesh are 
added to points where concentration gradient of the impurity material in 
the post-deformation diffusion calculation triangular mesh is large (step 
S80). 
Subsequently, after the setting of impurity material concentrations after 
deformation, another diffusion calculation triangular mesh for deformed 
shapes of regions of the semiconductor device is newly generated by the 
process simulator using the vertexes which has been added in the step S80, 
guaranteeing the Delaunay partitioning (new diffusion calculation 
triangular mesh) (step S81). 
Subsequently, new diffusion calculation control volumes corresponding to 
each vertex of the new diffusion calculation triangular mesh are defined 
by the process simulator (step S82). 
Subsequently, the amount of the impurity material in the post-deformation 
diffusion calculation control volumes is fully transferred to the new 
diffusion calculation control volumes by the process simulator by 
executing AND operation between the post-deformation diffusion calculation 
control volumes and the new diffusion calculation control volumes (step 
S83). 
Subsequently, diffusion calculation is executed by the process simulator 
using concentrations of the impurity material on vertexes of the new 
diffusion calculation triangular mesh and the new diffusion calculation 
control volumes, i.e. diffusion equation of the impurity material is 
solved (step S84). 
Subsequently, it is judged whether or not predetermined time has elapsed. 
If NO, process is returned to the step S74 to repeat the step S74 through 
the step S84 until the predetermined time is counted by the process 
simulator. Incidentally, in the next loop after the process is returned to 
the step S74, the new diffusion calculation triangular mesh in the 
previous loop is used as pre-deformation diffusion calculation triangular 
mesh in the step S78 of the next loop. If YES, the process simulation is 
ended (step S85). 
In the addition of vertexes of the oxidation calculation triangular mesh in 
the step S75, stresses of two adjacent vertexes of the oxidation 
calculation triangular mesh are compared and a vertex is added if the 
difference or the ratio between the two stresses is larger than a 
predetermined value. For example, it is possible to add vertexes to the 
midpoints of adjacent pairs of vertexes of the oxidation calculation 
triangular mesh if the stress ratio between the pair is larger than two. 
In the addition of vertexes of the diffusion calculation triangular mesh in 
the step S80, impurity material concentrations of two adjacent vertexes of 
the diffusion calculation triangular mesh are compared and a vertex is 
added if the difference or the ratio between the two impurity material 
concentrations is larger than a predetermined value. 
For example, it is possible to add vertexes to the midpoints of adjacent 
pairs of vertexes of the diffusion calculation triangular mesh if the 
ratio between impurity material concentrations of the pair is larger than 
ten. 
Incidentally, although vertexes are added in the step S75 and the step S80, 
it is also possible to reduce the number of vertexes on the oxidation 
calculation triangular mesh or the diffusion calculation triangular mesh 
in order to shorten calculation time. For example, it is possible to 
delete vertexes on the oxidation calculation triangular mesh where stress 
gradient is lower than predetermined value, or delete vertexes on the 
diffusion calculation triangular mesh where concentration gradient of the 
impurity material is lower than predetermined value. 
As described above, according to the second embodiment, it is possible to 
improve precision of oxidation calculation and diffusion calculation by 
adding vertexes of the oxidation calculation triangular mesh to points 
where stress gradient is large and by adding vertexes of the diffusion 
calculation triangular mesh to points where concentration gradient of the 
impurity material is large. 
As set forth hereinabove, by the process simulation method according to the 
present invention, use of two triangular mesh, i.e. the oxidation 
calculation triangular mesh which is suitable for oxidation calculation 
and the diffusion calculation triangular mesh which is suitable for 
diffusion calculation, is made possible, and process simulation keeping on 
guaranteeing the Delaunay partitioning and keeping conservation of the 
gross dose of the impurity material due to introduction such as ion 
implantation is realized even if shapes of regions of the semiconductor 
device vary as time passes due to oxidation, with both high precision 
calculation and short calculation time. Calculation time for the process 
simulation is reduced to half in comparison with the case where single 
triangular mesh is used. 
Furthermore, by adding vertexes of the oxidation calculation triangular 
mesh to points where stress gradient is large or by adding vertexes of the 
diffusion calculation triangular mesh to points where concentration 
gradient of the impurity material is large, precision of oxidation 
calculation or precision of diffusion calculation can be improved. 
While the present invention has been described with reference to the 
particular illustrative embodiments, it is not to be restricted by those 
embodiments but only by the appended claims. It is to be appreciated that 
those skilled in the art can change or modify the embodiments without 
departing from the scope and spirit of the present invention.