JPH09320919A - Manufacture of semiconductor device, manufacturing device, simulation method and simulator - Google Patents

Abstract

PROBLEM TO BE SOLVED: To provide a manufacturing method of a semiconductor device which enables a semiconductor device manufacturing process to be proceeded without a test piece as a desired process or while amendments being effected. SOLUTION: In a manufacturing method of a semiconductor device consisting of a plurality of processes, actual observation data in at least one of the plurality of processes is obtained. Prediction data in at least one of the plurality of processes is obtained by an abinitio molecular dynamics process simulator or a molecular dynamics simulator provided with empirical potential. Then, the prediction data and the actual observation data are compared and tested successively in an actual time. When a significant difference is recognized between a set value of a manufacturing process factor and a factor of the plurality of manufacturing processes gathered from the actual observation data according to the comparison and test, the manufacturing process factor is amended successively in an actual time.

Description

Detailed Description of the Invention

[0001]

BACKGROUND OF THE INVENTION 1. Field of the Invention The present invention relates to a semiconductor device manufacturing method, a manufacturing apparatus, a simulation method, and a simulator, and more particularly, to a limited quantity or qualitative content in a single-wafer type or clustered semiconductor device manufacturing apparatus. The method for manufacturing a semiconductor device and the method for manufacturing the semiconductor device, which maximizes the “in-situ” measurement value of the device and enables the semiconductor device manufacturing process to proceed according to a desired process or while being corrected without using a test piece The present invention relates to a device, a simulation method, and a simulator.

[0002]

2. Description of the Related Art Conventionally, semiconductor device manufacturing apparatuses have come to the present after undergoing the following changes. Looking at this in the flow of each generation, the dRAM generation and the chip cost for each generation change as shown in FIG. 128 as the capacity increases from 64K to 256M. FIG.
The vertical axis of 28 represents the chip cost at 1M as 1. Further, in FIG. 128, the plotted points show the chip cost, and the horizontal axis shows the dRAM generation.

As shown in FIG. 128, from 64K to 256
Up to K, the number of processes is relatively stable, and the chip cost is also stable. However, the number of processes is increasing every generation from 4M, and the chip cost is also increasing accordingly. In order to reduce the increase in cost due to the increase in the number of steps as much as possible, as shown in FIG. 128, for example, in 64M, (1) elimination of waste, (2) standardization of work, and (3) chip size Reduction, (4) Promotion of miniaturization,
Costs are being reduced by taking measures such as (5) productivity improvement, (6) FA promotion, and (7) larger diameter.

[0004] At 256M, by applying a new material, adopting a high dielectric material, applying a new structure, three-dimensional, etc., it will be put on the birule. However, when it comes to 1G class, the situation changes and it becomes unclear whether or not to follow the bi-rule. It is no longer possible to find a solution as a business by the conventional method that pursues the maximum performance of individual manufacturing equipment up to the 256M class.
Pursuing the highest performance of each manufacturing apparatus to the utmost means increasing the cost.

Further, as shown in FIG. 129, a device having high performance is provided for each process, for example, processes A to B → C.
It is assumed that they are used in the order of → D → E → F. In this case, it is obvious that the increase of the equipment area also causes the cost increase.

Another thing that must be forgotten in the 1G class and above is that precise control due to physical necessity becomes essential in each process. For example, it is the prevention of contamination, the elimination of natural oxide film, and the pursuit of planarization at the atomic level. At present, a clustered manufacturing apparatus as shown in FIG. 130 has been proposed as a new concept in consideration of these things, although it is still in the research stage. This clustered manufacturing apparatus (cluster tool) is of a batch type, but is mainly of a single-wafer type. Further, as can be seen from FIG. 130, for example, there is a transport system in the center, and each process is
As shown by the arrow, the sequence proceeds from A to B → C → D, respectively. At this time, in each step and in each step, in many cases, high vacuum is maintained, deposition of a natural oxide film can be reduced, contamination can be blocked, and surface control can be performed at the atomic level.

On the other hand, a field effect MIS formed on a semiconductor
With respect to the element, in the recent technological development of high speed and high integration, the ferroelectric is moving toward thin film. When the ferroelectric film is thinned, the threshold voltage becomes shallow, of course, and the main operation speed becomes faster by this amount, and thereby the AC characteristic is particularly improved. Also, considering an EEPROM or the like, as the device is miniaturized, it will be in a very severe usage condition. In such a case, the conventional ferroelectric film produced by a conventional method cannot provide sufficient reliability.

However, even if the MIS element is miniaturized, measures for improving the characteristics of the ferroelectric film itself are not neglected so much, and in the state where the power supply voltage does not drop so much, particularly in the main operation, the gate ferroelectric film is formed. A high electric field is applied. Further, electron holes generated by impact ionization or the like from the channel region are injected into the ferroelectric film depending on the boundary conditions such as the polarity of the gate electrode and the drain voltage. Then, these carriers are trapped in the ferroelectric film and not only deteriorate long-term reliability, but also cause decrease in breakdown voltage.

Let's look a little more at the atomic level. In the thermal ferroelectric film, for example, when a high electric field is applied to the silicon ferroelectric film, the Si-O bond forming the network becomes
It interacts with a high electric field applied from the outside. Then, the bond is broken and a trap center for trapping electrons and holes is formed. Electrons and holes that pass subsequently are trapped at the trap centers, and the electric field strength distribution in the film thickness direction is locally localized. It becomes higher than the average electric field, and eventually dielectric breakdown occurs.

[0010] In an attempt to improve these situations, the idea of single-crystallizing a ferroelectric film has recently been proposed at the academic level. For example, silicon (111)
Single crystal growth of cerium oxide CeO ₂ on the surface is reported in J. Appl. Phys., Aol. 69 (12), p8313 (1991). Alternatively, the fact that calcium fluoride CaF2 grows as a single crystal on a silicon single crystal is described in Japan.J.Appl.Phys.
Suppl., Aol. 21-1, p187 (1982). These are only a part, and not only are they still in the state of idea, but there are still many questions about the calculation method that is the basis of the guideline. Therefore, for example, even if a ferroelectric film is taken as an example, the single crystal of the ferroelectric film using it is not fully recognized. Let's look at what we have reported about the structure of this single crystal ferroelectric film. One is A. Miyamot, K. Takeichi, T. Hattori, M. Kubo, and
T. Inui, "Mechanism of Layer-Layer Homoepitaxial
Growth of SrTiO3 (100) as Investigated by Molecular
Dynamics and Computer Graphics ". Jpn. J. Appl. Ph
ys. Vol 31, (1992) 4463-4464. Here, the perovskite is used as a single crystal ferroelectric film, and it is used as a base S
The stable structure is calculated when it is bonded to the i-substrate. However, the initial arrangement of the angle or distance between Si—O—Ti or Ti—Si—Sr is erroneously used. Further, these simplify the structure of the single crystal ferroelectric film, and the guidelines are not always correct.

As can be seen from these facts, the present situation is that most of the discussion has not been made about how to design a single crystal ferroelectric film. On the other hand, with respect to the ferroelectric film itself, the present situation is that the improvement of the characteristics is still being refined in the forming process. For example, prepare a surface that is as clean as possible,
Alternatively, the discussion is all about the sputtering temperature and the sputtering atmosphere for film formation. It must be said that recognition of these, including single crystal ferroelectric films and structures related to ferroelectric films, is still extremely low. Moreover, the reason for the spontaneous polarization of the ferroelectric substance and its relation to the crystal structure are not known, so it is only a phenomenological pursuit.

[0012]

However, in the above-mentioned new cluster tools, the greatest difficulty currently facing is the situation observation method. That is, in the case of an apparatus that performs a single process, a wafer called a test piece was first introduced into each process. Using this test piece, for example, the film thickness, the refractive index, etc. were measured to test whether or not the process proceeded according to predetermined conditions.

Monitoring the conditions in each cluster tool by using such a test piece has a problem in quantitative and qualitative understanding. For example, an XPS or FT-IR apparatus is usually used in the oxide film forming step or growing step. Even if this is measured in-situ, the meaning of the data has not been sufficiently examined at present. At most, it is judged from the peak wave number position of the spectrum whether or not a predetermined substance is generated. In this situation,
Many rely on the judgment of researchers. In addition, it is not possible to respond quantitatively, and it is also not possible to respond promptly. Also,
When there is a new unknown substance that cannot be relied on by the experience of researchers, there is no specification for judgment.

On the other hand, the use of a molecular dynamics simulator has also been proposed. However, in general, the molecular dynamics simulator is (1)
It is said that the calculation is very slow, and (2) the target material execution area is very small. Therefore, it was almost unthinkable to directly connect the molecular dynamics simulator and the process variation calculation program. That is, since the molecular dynamics part is extremely slow to calculate,
Even if they were put together, it was thought that the whole thing would eventually be pulled to the slow speed of the molecular dynamics part, and would not be suitable for practical use.
In the molecular dynamics part, the calculation area is extremely small,
On the other hand, it is about several tens of angstroms at most. On the other hand, since the process variation calculation program calculates the area of several microns, it is difficult to directly connect them in terms of the deviation of the calculation area.

For example, the present inventors pay attention to the free energy of the entire system including the interface as the installation condition of the single crystal ferroelectric film, and position it so that this energy is minimized under the applied condition. MIS of high reliability
The device was made and the level was not formed at the interface. In particular, in determining the energy of the entire system, it is necessary to reliably reproduce the structure of the single crystal ferroelectric film and to grasp exactly how each of the constituent atoms moves. However, the present inventors have found that ab-ini
We completed tio molecular dynamics and established a unique potential integration method in molecular dynamics, and made these calculations possible for the first time. Moreover, it is intended to clarify the degree of quality of the single crystal ferroelectric film from the viewpoint of sufficiently assuring the improvement of the electrical characteristics and the like by using this.

For example, when a high electric field is applied to the tunnel insulating film, the tunnel insulating film interacts with a bond (for example, a Si-O bond in the case of a silicon ferroelectric film) constituting the basic structure of the tunnel insulating film, The bond is broken, and the broken bond is followed by the trapping of electrons and holes passing through the tunnel insulating film, and the electric field locally becomes larger than the average electric field depending on the charge distribution in the film and deteriorates. Progresses further. To suppress the progress of deterioration, the present inventor suppresses the breakage of the bond (for example, Si—O bond) of the insulating film due to the electric field.

The interaction between the electric field and the coupling of the insulating film depends on the direction of the electric field and the coupling direction of the insulating film. Therefore,
In order to weaken the interaction, it is necessary to reduce the number of insulating film bonding directions that are strong. When it is made completely random, such as amorphous, some of the bonds of the insulating film always have a strong interaction with the electric field.

An object of the present invention is to control the coupling direction of the ferroelectric insulating film, suppress the interaction with the electric field, and suppress the deterioration of the dielectric property due to the electric field.

This relates to a technique for producing a highly reliable dielectric film, which was created by the present inventors based on the calculation results by the new ab initio molecular dynamics system and original experimental data. If this is applied to a ferroelectric film, etc., the single crystal ferroelectric film is
If the positional relationship between the substrate and the substrate electrode is not optimized, it may be worse than the amorphous dielectric film. The criteria for the material and crystal orientation at this time are shown. As a new guideline by the inventor of the present invention, a simple word is "to arrange the crystal so that the free energy of the system is minimized under the main applied electric field". In addition, as a practical problem, "allowable range of crystal defect concentration that can exist in single crystal ferroelectric film" and the like are also described. In addition to the arrangement in which the free energy of the system at the time of main application becomes the minimum, the allowable range of unevenness of the free energy distribution when viewed locally in order to suppress the trigger of deterioration is also described.

Even if one single-crystal ferroelectric film is taken up, it has a very complicated structure such as a four-fold symmetric perovskite, and the characteristics are greatly changed depending on the combination of the bonding surface with the substrate. Although complicated, the inventors of the present invention try to "gain" under the main operating voltage even if "loss" occurs crystallographically.

Further, the present invention solves the problem of the monitor function in clustering, and enables "in-situ" measurement of a limited quantity or qualitative content in a semiconductor device manufacturing apparatus, particularly a single wafer type or clustering manufacturing apparatus. Provided are a semiconductor device manufacturing method, a manufacturing device, a simulation method, and a simulator that allow the semiconductor device manufacturing process to be expanded as much as possible without a test piece and according to a desired process or while being corrected. The purpose is to do.

[0022]

The present inventors have found that abin
By using itio molecular dynamics and developing a new variational calculation method, (i) the fluctuation of each process element and the superposition phenomenon of the deviation of each element are diagnosed, and (ii) whether or not the desired sequence is executed. We have developed a simulation system that enables the checking of (iii) and the discovery of (iii) fluctuation factors. Further features of this system are: (vi) process fluctuations can be handled mathematically and statistically and efficiently; (v) a new inference engine is added;
i) A new backward system is added.

That is, according to the present invention (claim 1), in a method for manufacturing a semiconductor device comprising a plurality of steps, a film formed in a manufacturing apparatus in at least one of the plurality of steps is formed on a semiconductor substrate. A film, a surface of a semiconductor substrate, and at least one of a plurality of properties of the inside of the semiconductor substrate to obtain actual observation data, and an ab initio molecular dynamics process simulator or empirical potential-added molecular dynamics By a science simulator, in at least one of the plurality of steps,
At least one set value of a plurality of manufacturing process factors, and a fluctuation of the set value as an input value, performing variational calculation based on the input value, and obtaining at least one predicted data of the plurality of characteristics; A step of sequentially comparing and testing the predicted data and the actual observation data in real time; at least one set value of the plurality of manufacturing process factors by the comparison test; and a plurality of the plurality of pieces estimated from the actual observation data. And a step of sequentially correcting the manufacturing process factor in real time when a significant difference is found between the manufacturing process factor and at least one manufacturing process factor.

According to the present invention (claim 2), in the method of manufacturing a semiconductor device (claim 1), the abiniti is used.
o A molecular dynamics process simulator or a molecular dynamics simulator given an empirical potential is a semiconductor material, an insulating material, or a semiconductor material formed or processed in the plurality of steps.
Alternatively, when the motion of the atoms forming the conductive material reaches equilibrium, the time is set to 0, and then the sampling time is set until the time t, and the polarizability, charge, and position vector of each atom are correlated with time. The prediction data derived from the ab initio molecular dynamics process simulator or the molecular dynamics simulator to which an empirical potential is given has a function of obtaining and obtaining an optical spectrum at a predetermined time. Provided is a method for manufacturing a semiconductor device, which includes the vibration behavior to be generated, its natural frequency, and distribution information of the natural frequency.

According to the present invention (claim 3), in the method for manufacturing a semiconductor device (claim 1), the abiniti
o A molecular dynamics process simulator or a molecular dynamics simulator given an empirical potential is a semiconductor material, an insulating material, or a semiconductor material formed or processed in the plurality of steps.
Alternatively, it has a function of sequentially calculating the material strength constant of the material from the motion of atoms constituting the conductive material, and using this to find the stress induced in the material, and when the stress exceeds a specified value. Provides a method for manufacturing a semiconductor device, wherein at least one of the manufacturing process factors is sequentially corrected.

According to the present invention (claim 4), in the method for manufacturing a semiconductor device (claims 2 and 3), the atoms constituting the semiconductor material, insulating material or conductive material are silicon atoms and oxygen atoms. A method for manufacturing a semiconductor device is provided.

According to a fifth aspect of the present invention, in a semiconductor device manufacturing apparatus for performing a plurality of steps, a film formed in the manufacturing apparatus in at least one of the plurality of steps is formed on a semiconductor substrate. Means for measuring at least one of a plurality of characteristics relating to at least one of the film, the surface of the semiconductor substrate, and the inside of the semiconductor substrate to obtain actual observation data;
An initio molecular dynamics process simulator or a molecular dynamics simulator given an empirical potential is used to input at least one set value of a plurality of manufacturing process factors in at least one of the plurality of steps, and a fluctuation of this set value as an input value. And, means for obtaining at least one prediction data of the plurality of characteristics by performing variational calculation based on this input value, the prediction data and the actual observation data are sequentially
A means for performing a comparative test in real time, and a significant difference between at least one set value of the plurality of manufacturing process factors and at least one of the plurality of manufacturing process factors estimated from the actual observation data by the comparative test. Provided is a device for manufacturing a semiconductor device, comprising means for sequentially correcting the manufacturing process factor in real time when a difference is recognized.

The present invention (claim 6) provides an abinitio
At least one set value of a plurality of manufacturing process factors in at least one of a plurality of steps of a manufacturing process of a semiconductor device by a molecular dynamics process simulator or a molecular dynamics simulator given an empirical potential, and Fluctuation is used as an input value, variation calculation is performed based on this input value, a film formed in the manufacturing apparatus in at least one of the plurality of steps, a film formed on the semiconductor substrate, a semiconductor substrate surface, and Obtaining at least one prediction data of a plurality of characteristics relating to at least one inside the semiconductor substrate, and at least one of the plurality of characteristics in the prediction data and at least one of the plurality of steps.
The present invention provides a simulation method, which comprises a step of diagnosing and judging the manufacturing process factor by sequentially performing real-time comparison and verification with real observation data obtained by measuring one of the two (claim 7). Is at least one set value of the plurality of manufacturing process factors by the comparison test and at least one of the plurality of manufacturing process factors estimated from the actual observation data in the simulation method (claim 6).
There is provided a simulation method comprising a step of obtaining a corrected value of the manufacturing process factor in real time when a significant difference is recognized between the two.

The present invention (Claim 8) provides a plurality of characteristics relating to the inside of the manufacturing apparatus, the surface of the semiconductor substrate, and / or the film formed on the semiconductor substrate in at least one of the plurality of steps of the manufacturing process of the semiconductor device. By using at least one of the plurality of manufacturing process factors in at least one of the plurality of processes, a means for measuring at least one to obtain actual observation data and an ab initio molecular dynamics process simulator or an empirical potential-provided molecular dynamics simulator. At least one set value and a fluctuation of this set value are used as input values, variational calculation is performed based on these input values, means for obtaining at least one predicted data of the plurality of characteristics, the predicted data and actual observation And a means for sequentially comparing and verifying the data in real time to diagnose the manufacturing process factor of the semiconductor device.
Provide a judgment simulator.

The present invention (Claim 9) is based on abinitio
At least one set value of a plurality of manufacturing process factors in at least one of a plurality of steps of a manufacturing process of a semiconductor device by a molecular dynamics process simulator or a molecular dynamics simulator given an empirical potential, and Fluctuation is used as an input value, variation calculation is performed based on this input value, a film formed in the manufacturing apparatus in at least one of the plurality of steps, a film formed on the semiconductor substrate, a semiconductor substrate surface, and Means for obtaining at least one predictive data of a plurality of characteristics relating to at least one inside the semiconductor substrate, and the predictive data and at least one of the plurality of characteristics in at least one of the plurality of steps.
And sequentially observing the real observation data obtained by measuring one of them in real time, and diagnosing the manufacturing process factor of the semiconductor device, comprising:
Provide a judgment simulator.

On the other hand, the inventors of the present invention have proposed new ideas for the molecular dynamics simulator with respect to the problem that the calculation is very slow and the problem that the execution area of the target substance is very small. It was introduced and overcome, and the molecular dynamics simulator and the process variation calculation program were combined by a new method, and by this combination, it was possible to bring out a new function that has never existed before.
This will be described below in order.

First, the problem of calculation speed (1) will be touched on in some detail. In molecular dynamics, a motion equation is set up for each atom according to, so to speak, Newtonian mechanics, and the equations of motion are solved simultaneously over all 100 or 1000 atoms of interest.

The calculation of the force that is the basis of this equation of motion is the most complicated. This is because the conventional calculation method is slow. Each atom receives force from the atoms around it. It receives power from both nearby and distant atoms. This needs to be taken into consideration. Generally, the potential energy V _ij derived according to the Telesoff formula is obtained, and the force is calculated based on this potential energy V _ij . However, in the past, in general, a method of preparing a table for each direction of the force applied to each atom in advance based on the potential, the distance, the direction, and the like was generally used. In this method, the magnitude of the potential is checked for each 100 or 1000 atoms by referring to the table every time. This is very time consuming. In addition, the force cannot be calculated from the potential itself, and the inclination of the magnitude of the potential with respect to each direction is
It consists of forces in each direction (see Fig. 2). Therefore, it can be seen that it takes a lot of computer time.

The reason why the numerical table is created in this way is that the equation representing the shape of the potential is very complicated. As mentioned above, potential (V)
The positions of the surrounding atoms are complicatedly entangled with each other in the angle (θ) and further in the distance (r). Therefore, it is very difficult to differentiate it mathematically to create a function. For example, θ is a function of r, but it turns around and becomes a function of r again.

The inventors of the present invention are very complicated in general as described above, and in the case where the argument of the potential expression is cyclic as described above, in the part where the partial differentiation has not been attempted in the past, a new A higher-order cyclic partial differential equation was created. This makes the transformation of the equation itself seemingly complicated, but the calculation can be done very quickly. As before,
I no longer need to refer to the number table.

Then, the first concept of this higher-order cyclic partial differential equation is shown below. This generally assumes 3 atoms,
thinking. That is, if three atoms are generally described, then this method can be expanded. If there are three atoms, the distance to each atom and the angle between the atoms of interest can also be calculated. From these, the potential is uniquely obtained, and based on this, the “slope” of the potential can be easily calculated in any direction. The present inventors have confirmed that the differential type is correctly created by creating partial differential equations for all the variables that are cyclically used and solving them in parallel.

The outline of the system according to the present invention will be specifically described below. First, the actual application of the present invention will be described with reference to FIG. That is, the process device V is a part of the semiconductor manufacturing system which is actually operating. here,
The absorption spectrum of light is observed by the spectrum meter VI.

The light absorption spectrum obtained by the spectrum meter VI is once stored as a file 4,
Sent to scheduler VII. Predictor values according to the present invention, which will be described later, are also input to the scheduler VII, and these are compared by the comparison diagnostic test unit IV.

If the comparison result in the comparison diagnosis verification unit IV exceeds the allowable value, it is determined that the specification is not satisfied and the line is removed from the line. If satisfied, the variation is sent to the process variation calculation unit I as an inference engine.

The process variation calculation unit 1 basically calculates a macroscopic property not considering the atomic level and stores it in the file 2. Regarding the contents of file 2, the interatomic position vector and bond angle are calculated by the molecular dynamics simulator II.

The interatomic position vector and bond angle calculated by the molecular dynamics simulator II are obtained by the film property calculation unit I.
II, where the basic film properties are calculated and returned to File 3. The basic film properties calculated by the film property calculation unit III are also returned to the process variation calculation unit I to be a part of the parameters.

As described above, according to the present invention, a method of manufacturing a semiconductor device is realized which allows a semiconductor device manufacturing process to proceed without a test piece, according to a desired process or while being corrected.

In practice, the content of the molecular dynamics simulator II is particularly important. Specifically, as shown in the block diagram of FIG. 7, it is composed of two parts. The first part IIa uses the higher-order cyclic partial differential method introduced by the inventor of the present invention, and the second part IIb performs analysis considering higher-order terms of third or higher order. The details of the higher-order cyclic partial differential method and the analysis considering higher-order terms of the third or higher order will be sequentially described below.

That is, the present inventors have returned to pure mathematical theory, established a new variational calculation method, and developed a technique for predicting various spectra by a new computational chemistry theory. Hereinafter, the new abiniti in the present invention
o Explain molecular kinematics.

The overall structure of the new simulation developed by the present inventors is shown in FIG. Molecular dynamics part (MD)
We have proposed a new method that synthesizes and the density functional (DFT). First, the use of the molecular dynamics method and the density functional theory method by the present inventors will be shown below with reference to FIG.

As shown in FIG. 3, the molecular dynamics (MD) part occupies the lower half of the flowchart, and the density functional part (DFT) occupies the upper half. Molecular dynamics method is zero K
It deals with a rather large system at a so-called finite temperature, and uses the potential obtained in advance from the first principle as the potential between atoms (ions). On the other hand, density functional (DFT)
Was used to find the ground state of a small system.
And, naturally, the temperature of the DFT portion was set to 0K. The potential part was obtained from the ground state of the electronic system. Specifically, the wave function (KS equation) was solved and the local density functional was used together.

At the DFT portion, the total energy of the ground state of the electronic system without degeneracy is determined and set. The total energy is the functional E of the electron density n (r) (without going back to the wave function)
[N (r)].

That is, E [n (r)] is expressed by the following equation (120)
It is represented by 1).

[0049]

[Equation 1] Vz (r) of the first term on the right side is an external field such as Coulomb field of the nucleus,
The second term is the Coulomb interaction energy between the electrons, and the third term is the kinetic energy Ts [n (r)] when there is no electron interaction, assuming the one-body approximation of the electrons. 1
202).

[0050]

[Equation 2] The fourth term is the exchange correlation energy Exc [n (r)], the shape of which depends on the local density approximation. The electron density n (r) is represented by the following formula (1203).

[0051]

(Equation 3) As the ground state, KS shown in the following equations (21) and (22) can be obtained by setting the variation of E [n (r)] of the equation (1201) as 0.
The equation (Kohn-sham equation) is solved, and N number of eigenvalues E is adopted. Expression (12
Since n (r) is included in 11), n (r), ψ
(R), the effective potential Veff (r) must be solved self-consistently.

[0052]

[−∇ ² + V _eff (r)] Ψ (r) = E Ψ (r) (1211) V _eff (r) = V _Z (r) + ∫2n (r ′) / (| r− r ′ |) · dr ′ + δE _xc [n (r)] / δn (r) (1212) The exchange correlation energy Exc [n (r)] is the energy per electron εxc (in the local density approximation (LDA). n) is multiplied by the electron density n and integrated to obtain the following formula (121
3) was obtained.

[0053]

[Equation 5] E _xc [n (r)] = ∫ε _xc (n (r)) n (r) dr (1213) Various forms of εxc (n) have been proposed, such as Wigner et al. . Furthermore, when treating only valence electrons without considering inner-shell electrons, the Coulomb field of the nucleus is not used as it is as Vz (r) in equation (1211), but a norm-conserving pseudopotential that eliminates the singularity at the nuclear position is used. Change. By smoothing in this way, the number of components in the Fourier expansion could be suppressed.

As the total energy, the present inventors have used the one including Coulomb potential between atoms (ions) for the sake of strictness. That is, it is the following formula (1221).

(Equation 6) Where {R _I} nuclear coordinates of ions, {αν} external constraints (volume Omega, etc. strain εLν). Formula (122
The first term on the right side of 1) is Ts [n (r)] on the right side of the DFT equation (1201), and the rest of the right side of equation (1201) is a part of U of equation (1221). That part is DF
Similar to the case of T, V _eff
It is represented by (r), and Vz (r) is replaced with a pseudopotential and used.

A new way of thinking is entered here. That is, E is regarded as potential energy, and in addition, virtual kinetic energy shown in the following equation (1222) is introduced.

[0056]

(Equation 7) Therefore, the Lagrangian is represented by the following formula (1223).

L = K−E (1223) In formula (1222), M _I is the mass of a real atom (ion), but μ and μν are virtual masses, and are, as described later, depending on the purpose. changed. As a constraint condition, the following formula (122)
Orthonormal condition shown in 4)

(8) Since ∫ψ _i ^* (r, t) ψ _k (r, t) dr = δ _ik (1224) is imposed, an undetermined multiplier Λ _{ik is generated} when the Euler equation is _created.
Is introduced, and the variation of the sum of L and the formula (1231) below is taken. As a result, the following equations (1232) and (123
3) and (1234) are obtained.

[0058]

[Equation 9] The temperature T of the kinetic energy in equation (1222) corresponds.
Because of the degree of freedom by the energy equipartition law of classical statistical mechanics have temperature, in terms of the relationship between _{Σ1 / 2 M I R I 2} , the temperature of the real rather than virtual. Also μ, μ
ψi and αν can be suppressed or enhanced depending on the size of ν.

For example, if μ << M _I and T → 0 from a high temperature, ψ i changes while R _I is stopped and the ground state of the electronic system is obtained. At this time, the left side of the equation (1231) becomes 0 and the right side becomes the following equation (30) due to the equation (1221) and the cautions thereunder. Therefore, {Λ _ik } that appropriately forms a linear combination of ψ i is a diagonal line. To obtain the KS equation (1235).

## EQU10 ## [∇ ² −V _eff (r)] Ψ _j (r) + ΣΨ _k (r) (1235) In other words, the same result as solving the KS equation with DFT is obtained.
It was obtained through the muted annealing. However, since the state of the temperature T is created by following the motion of the virtual dynamical system regardless of the Monte Carlo method, d
dynamically simulated anneali
ng (DSA).

When the equation (1231) is integrated, ψi is not necessarily treated as a variable as it is, but the expansion coefficient (that is, Fourier coefficient) obtained by expanding it with a plane wave can be treated as a variable. In order to prevent the number of expansion coefficients from unnecessarily increasing, it is necessary to eliminate the singularity at the nuclear position by pseudopotential as in the case of DFT. In particular, when the periodicity of the crystal can be utilized, the linear combination of ψ i that makes {Λ _ik } diagonal is described as φnk.
k is a wave number vector in the Brillouin band, and n is a band index.

Φnk should satisfy the following expression (1241).

[0062]

[Mathematical formula-see original document] μS _nk (r, t) = [∇ ² −V _eff (r)] φ _nk (r) + λ _nk φ _nk (r, t) (1241) λnk is an eigenvalue of the matrix {Λ _ik } Yes, it is an energy level. φnk (r, t) is expanded and the following formula (1242)
Is obtained.

[0063]

(Equation 12) If the pseudopotential is local (this condition is not generally satisfied), the effective potential V _eff is also expanded to obtain the following formula (1243).

[0064]

(Equation 13) By substituting the equations (1242) and (1243) into the equation (1241), the following equation (1251) is obtained.

[0065]

[Equation 14] If K = K ′ + K ″ at the end of the right side, Σexp [i
Since (K + G) r] is common to all terms, the following equation (1252) is obtained.

[0066]

## EQU15 ## μT _{n, K + G} (t) = [− | K + G | ² + λ _nk ] C _{n, K + G} (t) −ΣV _GG ′ (t) C _{n, K + G} ′ (t) = ^{- [| K + G | 2} + V 0 (t) -λ nx] C n, K + G (t) -ΣV GG a _{'(t) C n, K} + G' (t) ... (1252) where the frequency, Define as shown in the following formula (1253),

Ω = ((| K + G | ² + V ₀ (t) −λ _nx ) / μ) ^1/2 (1253) When the left side of equation (1251) is corrected to the difference, the following equation (12)
61) is obtained.

[0067]

C _{n, K + G} (t + Δt) = 2cos (ωΔt) C _{n, K + G} (t) −C _{n, K + G} (t−Δt) −2 [1-cos (ωΔt)] { Σ _G ′ V _GG ′ (t) · C _{n, K + G} ′ (t)} (1261) Since this formula analytically integrates the oscillating part, it is better than a purely numerical method. The calculation amount can be reduced by increasing Δt.

The Coulomb potential is as follows. That is, the present inventor has found that the Coulomb potential is decomposed into four terms. In particular, it was confirmed that there was a new fourth term, although there were only three terms in the past. These methods will be sequentially described.

The inventors of the present invention first started to solve the strict formula from the basic formula in which the dielectric constant was taken into consideration. First, the solution was started from the following formula (1262), which is a basic formula.

[0070]

(Equation 18) The Coulomb potential is different between the conductor ε = ∞ and the vacuum ε = 1, L is one edge of the (cubic) unit crystal, Σ is in the unit crystal, and Z i and r i are the charges of the i-th particle, respectively. The position. This is because the charge inside the sphere causes polarization on the inner surface of the sphere. There is a layer of dipoles on the inner surface of the sphere that is not a conductor, but the last term in the above equation just works to counteract that effect. Ewald's method is on the left side, ε = ∞
Therefore, the last term of the above equation must be added to obtain the value in vacuum. Only the results are listed.

The Coulomb potential is expressed by the following equation (127)
Shall be represented by 1)

[Equation 19] In formula (1271), N is the number of atoms in the unit crystal, the position and charge of the first atom in them are ri and Zi, and n is a vector that specifies the unit crystal and its periodic shift. Therefore, it was set as shown in the following formula (1272).

It is expressed as n = n _x Z _x + n _y Z _y + n _z Z _z (1272). Where Zx, Zy, and Zz are x,
Vectors of edges in the y and z directions, and nx, ny, and nz are integers (in the case of bulk crystal) ranging from -∞ to + ∞. The symbol ′ in Σ ′ means that j = i is excluded when n = 0.

Here, the present inventors introduce a new F function shown in the following equation (1273).

[0074]

(Equation 20) In the above formula (1273), G (r, t) is a periodic function of t. We have found that this can be represented by a Fourier series.

G (r, t) can be further modified as shown in (1281) below.

(Equation 21) In the above equation (1281), m is a reciprocal lattice vector and is represented by the following equation (1282).

[0076]

[Number 22] _{m = m x · (1 /} L x, 0,0) + m y · (0,1 / L y, 0) + m z · (0,0,1 / L z) ... (1282) index mx, my and mx are integers ranging from -∞ to ∞, respectively. When m = 0, exp [...] = 1, so ΣZ
If i = 0 and the total charge in the unit crystal is 0, m =
The 0 term disappears. Excluding the term of m = 0 beforehand, the following equation (1283) is obtained.

[0077]

(Equation 23) The present inventors obtain the following equation (1291) by dividing the integration range by k and using two forms of G (r, t).

[0078]

(Equation 24) Here, the present inventors used the following formula (1292) as a formula.

[0079]

(Equation 25) As a result, the initial Coulomb potential is expressed by the following equation (1293).

[0080]

(Equation 26) The mark 'on the upper right of Σ'means that j = i is excluded when n = 0.

Further, the following equation (1301) is established.

[0082]

[Equation 27] Since this equation (1301) does not include r, it does not relate to force. This equation (1301) can be further transformed into the following equation (1302).

[0083]

[Equation 28] Summarizing the above results, the Coulomb potential is ultimately expressed by the following equations (1311) to (1315).

[0084]

(Equation 29) Here, n is n = nx * (Lx, 0,0) + ny *
(0, Ly, 0) + nz · (0, 0, Lz), and
m is m = mx * (1 / Lx, 0,0) + my * (0,
1 / Ly, 0) + mz * (0,0,1 / Lz).
The good point of the above expression is that the original Σ term attenuates only in the order of reciprocal, whereas the Σ term in Φ (1) depends on the factor of erfc, and the Σ term in Φ (2) The term is rapidly decayed by the factor of exp. Since it works in the opposite direction against the slowdown with Φ (3), it is necessary to select an appropriate κ that balances both. These are the results of calculating the contribution of the Coulomb force in order from the closer distance, and are the results when the periphery is a conductor.

Another term is added if the ambient is a vacuum. This is a part that has not been particularly omitted.

Generally, in molecular dynamics, the interatomic potential is the most important amount. At present, there are quite a few things that are simply described between two bodies. This is because the molecular dynamics calculation takes a very long time, and in order to save the time, a simple potential is used. In this way, if the potential is simplified, the calculation itself becomes faster, but it cannot reflect the actual situation.

As described above, the important interatomic potential can be generally described by three bodies. The implications of this are as follows. That is, there is an effect that the action of the first atom (i) and the second atom (j) enters from the third atom (k) through the second (j). This is shown in FIG. As shown in the figure, r is the distance between particles and θ is the particle angle. The master-slave relationship of these various quantities is as shown in FIG. In particular, it can be seen that it is very complex, as can be seen in FIG. Generally, the potential is differentiated to become a force, but as you can see, the calculation becomes very complicated. Here again, the conventional calculation method and the positioning of the present invention are shown using FIG.

The present inventors newly developed the higher-order cyclic partial differential method shown in FIG. 6 so that high-precision and high-speed operation can be performed. How fast the method of the present inventors is will be described below by including an actual comparative example. First, the simple two-body potential calculation shown in FIG. 6 will not be described here. This is because the description of the physical phenomenon of the semiconductor LSI cannot be obtained with sufficient accuracy to be used, so that it will not be discussed here. On the other hand, the present inventors have created a table reference method program as a conventional method with the potential of including the three-body problem. According to this, a table was created for r _ij , θ _{i jk} , and r _jk . Actually, it was stored as an array as block data. And, in consideration of the time when it deviated from this reference item, I also created a complementary program. The cubic spline method was used for complementation. The program is such that a table can be referred to for a certain value, that is, r _ij , θ _ijk , and r _jk . However, it is necessary to find out which matrix in the array the desired condition (r _ij , θ _ijk , and r _jk ) falls into. The present inventors
I wrote the program while paying attention to making it as fast as possible, but again, I needed a part to separate cases according to the size of the numbers, and I used an if then else statement there.
In this way, it was found out which frame of the table the desired condition (r _ij , θ _ijk , and r _jk ) fits in, and the complementation in it was also performed. With this, first, the potential under the desired conditions (r _ij , θ _ijk , and r _jk ) was obtained, but what is actually desired is the force at this point. The present inventors obtained the potential at the other point of the minute change using the same procedure, and from the difference between the two, obtained the derivative, that is, the force, using the idea of the average change rate.

As can be seen, there are so many if the
n else procedure was used. On computer, if then else
The procedure is relatively weak and slow. In fact, when the reference type program created by the present inventors was used, only the part for calculating the force was examined in detail. As a result, 160
Single crystal substrate consisting of 0 Si atoms, absolute temperature 300
In the case of K and an external pressure of 1 atm, the reference type program took about 350 times longer. As you can see,
It can be seen that the higher-order cyclic partial differential method proposed by the present inventors is more advantageous. In the reference program, it is necessary to create a data table for each atom.

Next, how to extract material characteristic values from molecular dynamics and how to extract them into a general-purpose simulator will be briefly described.

Usually, various evaluation techniques are used in the LSI process. Evaluation method of physical properties of the film by infrared absorption, further stress evaluation in the film by Raman spectroscopy,
There is a structural analysis of the film by X-rays. In addition to these optical techniques, there are also electrical evaluation techniques such as the dielectric constant of the film.

The inventors of the present invention have eagerly and theoretically decomposed the principle of each evaluation method into motions of individual atoms and traced them in detail. Then, by using the molecular dynamics developed by the present inventors, it is possible to trace the movement of individual atoms, theoretically calculate the evaluation amount of each atom for the first time, and directly display it in a form such as a spectrum. done. As a result, for example, a large amount of progress in technology that has never been achieved has been achieved, such as quantitative comparison with measured values.

The tracing of atomic motion by the theoretical means which the present inventors have overcome for the first time and the derivation of a concrete measurement evaluation amount will be described below. There are several, but basically,
The basic concept derived by the present inventor is as follows. That is, the moment of the atom of interest is calculated. Then, this is integrated within a window of time of interest, and then Fourier transformed to obtain an eigen spectrum. This is the main point.

Next, taking Si as an example, how the present inventor converted the potential into a force with respect to the empirical potential is shown. This is an example, but Si
If the effect of ionic coupling is not taken into consideration, the higher-order cyclic partial differential described here can be used as the three-body problem for the potential of the oxide film, for example, like the problem of Si. The present inventors also examined what kind of potential is good with respect to Si.

The potential used here is according to Tersoff. (Phys. Rev. Lett. 56.632 (1986), Phy
S. Rev. B37, 6991 (1988)).

Now, according to Tersoff, the total potential for the i-th Si is

[Equation 30] Can be described by Since this is described with three bodies, what the present inventors want to call attention to is V _ij ≠ V ji in the above equation (0001). If the position number of the Si atom of interest is i and other particle numbers around it are j, the above V _ij is

V _ij = f _c (r _ij ) [a _ij f _R (r _ij ) + b _ij f _A (r _ij )] (0002). Where r is the distance between the particles. Also, fc
(R _ij ) is a cutoff function, fR (r _ij ) shows repulsive force, and fA (r _ij ) shows attractive force. a _ij is a cutoff coefficient considering the coordination number, and b _ij is a cutoff coefficient considering the coordination number. 3 using the parameter expressing the coordination number
You are doing physical expression. Since the basic form is the same for other devices, the following description will be given here together with the idea of the traveling method developed by the present inventors for the first time.

FR (r _ij ) and fA (r _ij ) are modified Morse type potentials, and fR (r _ij ) = Aexp (-
λ1r) and fA (r _ij ) = − Bexp (−λ2r).

Of these, λ1 and λ2 are constants, and their magnitudes are the reciprocals of the values of the interatomic distance.

By the way, the cutoff function fc (r _ij ) is

Equation 32] _{V ij = f c (r ij} ) [a ij Aexp (-λ 1 r ij) -b ij Bexp (-λ 2 r ij)] f c (r) = 1 (r ≦ R-D) f _c (r) = 1 / 2−1 / 2 sin [π / 2 (r−R) / D] (R−D <r <R + D) f _c (r) = 0 (r ≧ R + D) (0003) Where R is sized so that it typically includes only the first adjacent zone of the structure of interest. The value is approximately 2 to 3 Å according to the detailed examination by the present inventor. Then b
_ij indicates the actual coordination number, but the cutoff function is used here as well. According to Tersoff, the definition is b _ij = (1 + β ⁿ ζ _ij ⁿ ) ^{−1 / 2n} (0004). here,

Ζ _ij = Σf _c (r _ik ) g (θ _ijk ) exp [λ ₃ ³ (r _ij −r _ik ) ³ ] ... (0005). The Σ symbol turns with k ≠ i, j. As you can see here,
Since the meaning of ζ _ij is an environmental factor due to the _inclusion of the third atom k, the size is different between i and j. That is, ζ _ij ≠ ζ jj, and further, V _ij ≠ Vji as described in the above equation (0001). Therefore, b _ij ≠ b ji.

Further, g (θ) is a bond angle factor,

G (θ) = 1− (c ² / d ² ) −c ² / (d ² + cos θ ² ) ... (0006) Here, θ is taken as shown in FIG. When finding θ, try to express it using actual Cartesian coordinates. That is,

Equation 35] r _ij = √ a ^{_{{(x j -x i) 2}} + (y j -y i) 2 + (z j -z i) 2} ... (0007), r ik in a similar procedure Desired. Then, if the inner product is P _ij ,

P _ijk = (x _j −x _i ) (x _k −x _i ) + (y _j −y _i ) (y _k −y _i ) + (z _j −z _i ) (z _k −z _i ) ... (0008). By using these, cos θ _ijk = P _ijk / (r _ij · r _ik ) (0009). Here, the constants in the above equations used by the present inventors are shown. That is,

Equation 37] R = 3.0Å, D = 0.2Å, A = 3264.7ev, B = 95.373ev, C = 4.8381, λ 1 = 3.2394Å, λ 2 = 1.3258Å, λ 3 = λ 2, β = 0.33675, n = 22.956, d = 2.0417 (0010).

After going through these preparations, the traveling part of the present inventor will be described in detail. From here, the power is finally calculated. Differentiating the potential (0002) with respect to the position coordinates gives a force. That is,

(38) Are the x components of the force vectors acting on the particles i and j, respectively. However, actually seeking it is not so easy. The inventions and inventions that return to the mathematical theory of the present inventors are as follows. That is the division of higher-order cyclic partial differentiation.

The expressions described in the equations (0011) and (0012) may be partially differentiated specifically. But it's easy to see that it's not that simple.

That is,

[Equation 39] Further, the modification of the partial differential equation with respect to j and k is also as follows. Further, although the above description mainly describes x, it is necessary to perform y and z. Here, in particular, the blank portion is left blank so that the correspondence with the above can be seen. The present inventors have not found such a complicated deformation as far as they have examined. Since we have never overcome such complicated transformations, we have instead escaped to a simple table reference method. By overcoming such an invention, the present inventors were able to combine molecular dynamics and a general-purpose fluid simulator.

[0104]

(Equation 40) Each term in these equations can be transformed as follows.

[0105]

[Expression 41] ∂V _ij / ∂r _ij = ∂f _c (r _ij ) / ∂r _ij・ [Aexp (−λ ₁ / r _ij ) −b _ij Bexp (−λ ₂ / r _ij )] + f _c (r _ij ) [− λ ₁ Aexp (−λ ₁ / r _ij ) + λ ₂ b _ij Bexp (-λ ₂ / r _ij )] = Aexp (−λ ₁ / r _ij ) · [∂f _c (r _ij ) _{/ ∂r ij -λ 2 f c (} r ij)] - Bexp (-λ 2 / r ij) · [∂f c (r ij) / ∂r ij -λ 2 f c (r ij)] b ij ... (0016) ∂f _c (r _ij ) / ∂r _ij = −π / 4D · cos [π / 2 · (r−R) / D] (R−D <r <R + D) ∂f _c (r _ij ) / ∂r _ij = 0 (r ≧ R + D) (0017) ∂V _ij / ∂ζ _ij = (∂V _ij / ∂b _ij ) ・ (∂b _ij / ∂ζ _ij ) = −Bfc (r _ij ) exp -(λ2r _ij ) (-1 / 2n) (1 + β ⁿ ζ _ij ⁿ ) − ^{1 / 2n-1} β ⁿ ζ _ij ^n-1 = Bfc (r _ij ) exp- (λ2r _ij ) b _ij (βζ _ij ) ⁿ / [2 {1+ (βζ _ij ) ⁿ } ζ _ij ] (0018) ∂ζ _ij / ∂r _ij = 3λ ₃ ³ Σf _c (r _ik ) g (θ _ijk ) exp (λ ₃ ³ (r _ij −r _ik) ³⁾ (r _ij - _ik) ² ... (0019)

[Number 42] _{∂ζ ij / ∂r ik = df c} (r ij) / dr ij · g (θ ijk) exp (λ 3 3 (r ij -r ik) 3) -3λ 3 3 f c (r ik ) f _c (r _ij ) (r _ij −r _ik ) ² … (0020) ∂ζ _ij / ∂cos θ _ijk = f _c (r _ik ) exp (λ ₃ ³ (r _ij −r _ik ) ³ ) dg (θ _{ijk / dcosθ ijk f c (r} ik) exp (λ 3 3 (r ij -r ik) 3) [2c 2 cosθ ijk / {d 2 cos 2 θ ijk} 2] ... (0021) ∂r ij / ∂x _i = (x _i −x _j ) / r _ij = ∂r _ij / ∂x _i ((0022) ∂r _ik / ∂x _i = (x _i −x _k ) / r _ik = ∂r _ij / ∂x _k … (0023)

∂cos θ _ijk / ∂x _i = 1 / (r _ij r _ik ) ・ ∂P _ijk / ∂x _i ⁺ P _ijk {1 / (r _ik ) ∂ / ∂x _i [1 / (r _ik ) ] ⁺ 1 / (r _ij ) ∂ / ∂x _i [1 / (r _ik )]} = 1 / (r _ij r _ik ) ・ (x _i −x _k + x _i −x _j ) −P _ijk {(x _i x _j ) / (r _ik r _ij ³ ) + (x _i −x _k ) / (r _ik r _ij ³ ) = 1 / (r _ik ) ・ ((x _i −x _j ) / (r _ij ) − (x _i −x _k ) / (r _ik ) cos θ _ij } ... (0024) ∂cos θ _ijk / ∂x _i = −1 / (r _ij ) ・ {(x _i −x _k ) / (r _ik ) − ( x _i −x _j ) / (r _ik ) cos θ _ij } ... (0025) ∂cos θ _ijk / ∂x _k = −1 / (r _ik ) ・ {(x _i −x _j ) / (r _ij ) − (x _i− x _k ) / (r _ik ) cos θ _ij } ... (0026) L (r _i ∂r _i / ∂t, V, ∂V / ∂t) = 1 / 2Σm (∂r _i / ∂t) ² − U (r _i ) +1/2 M (∂V / ∂t) ² −P _E V (∂L (q _j , q ′ _j ) / ∂q _j ) −d (∂L / ∂q ′ _j ) / dt = 0 ... (0027) FIGS. 7 and 8 show the molecular dynamics simulator II and The cooperation with the process variation calculation unit I will be described. In the molecular dynamics simulator II, the first input a is
The files (22), (32), (42), shown in FIG.
Some are transferred from (52), (62), (72), etc. The contents are boundary conditions, area size, stress, strain, temperature, and the like. If necessary, there is an auxiliary input b as the second one, which is used to specify the calculation time of molecular dynamics, the ensemble, and the like.

On the other hand, the position vector r (t) appears as the first output. Then, the files (23), (33), (43), (53), (63), shown in FIG.
(73), (83), etc. In addition, the film characteristic calculation program III calculates the film characteristics such as optical, electrical and crystallographic characteristics, and transfers them to the file (3) in FIG.

FIG. 9 shows an example of the input sequence of the general-purpose 2 / 3-dimensional process simulator ((Ia) in FIG. 8). Further, FIG. 10 shows another input sequence of the general-purpose 2 / 3-dimensional process simulator ((10) in FIG. 8). The contents shown in FIG. 9 also specify, for example, the respective cluster tools to be used, and the details of the oxidation process therein are, for example, the gas amount, the temperature raising sequence, and the valve sequence.

On the other hand, the input sequence of FIG. 10 is a case where the unknown material shown in FIG. 8 is included. In the case of FIG. 10, it is a ferroelectric film, for example, and includes a step of sputtering the film. At this time, since there is no detailed data of the ferroelectric film, jump to the molecular dynamics simulator II,
Calculation is performed there.

Further, the present inventors conducted a molecular dynamics simulation with an empirical potential, and (1) diagnosed fluctuations of each process element and superposition phenomenon of deviation of each element,
We have developed a simulation system that enables (2) checking whether the sequence is being executed according to the desired sequence and (3) discovering variation factors.

Specifically, the positional relationship between the single crystal ferroelectric film and the Si substrate and the positional relationship with the metal electrode are proposed, and further, the actual quality of the single crystal ferroelectric film is improved. For the first time, I was able to give a guideline as to whether or not to prepare something. The structure of such a single crystal ferroelectric film is reproduced by a computer using a rigorously unique calculation method, and a characteristic evaluation function is given. Based on this, the positional relationship with the substrate single crystal is considered in consideration of the application conditions. At the same time, the guideline for the allowable range of the defect density is shown.

That is, the present invention (claim 16) is a method of manufacturing a semiconductor device comprising a plurality of steps, wherein a film formed in a manufacturing apparatus in at least one of the plurality of steps is formed on a semiconductor substrate. Film, semiconductor substrate surface,
And a step of obtaining at least one of a plurality of characteristics relating to at least one of the inside of the semiconductor substrate to obtain actual observation data, and by a molecular dynamics simulator given an empirical potential, in at least one of the plurality of steps,
At least one set value of a plurality of manufacturing process factors, and the fluctuation of this setting as an input value, performing a variational calculation based on this input value, obtaining at least one prediction data of the plurality of characteristics, and Prediction data and actual observation data, a step of sequentially performing a comparison test in real time, and by the comparison test, at least one set value of the plurality of manufacturing process factors, and the plurality of inferred from the actual observation data A manufacturing method of a semiconductor device, comprising the step of sequentially correcting the manufacturing process factor in real time when a significant difference is recognized with at least one of the manufacturing process factor.

According to the present invention (claim 17), a molecular dynamics simulator provided with an empirical potential used in the method for manufacturing a semiconductor device described above is a semiconductor material, an insulating material, or a semiconductor material formed or processed in the plurality of steps. Alternatively, when the motions of the atoms forming the conductive material reach in parallel, the time is set to 0, and then the sampling time is up to the time t, and the time correlation with the charge and position vector of each atom is calculated, Predictive data derived from the molecular dynamics simulator given the empirical potential has a function of obtaining an optical spectrum in time, vibrational behavior derived from the motion of each atom at each time and time, its natural frequency,
Also provided is a method of manufacturing a semiconductor device including information on the distribution of natural frequencies.

According to the present invention (claim 18), the position vector of an atom and the electronegativity of each element, or each element, when calculating the motion of an atom in the molecular dynamics simulator given the above empirical potential, The method for manufacturing a semiconductor device is described in which the ionization potential and the electron affinity of are input, the change in charge due to vibration and displacement of individual atoms is calculated from these data, and the potential is calculated in consideration of this change in charge. provide.

By the new computational chemistry theory, the present inventors
We have invented a method for calculating and predicting various phenomena in a system containing various elements at high speed and with high accuracy. This has made it possible to calculate the properties of various materials used in the semiconductor manufacturing process with high accuracy.

Conventionally, as a means for predicting the structure and characteristics of crystals and molecules containing various elements, abinit
The io molecular dynamics method and the first principle molecular orbital method have been used. The ab initio molecular dynamics method is mainly used for investigating the structure and electronic state of a crystal, and the molecular orbital method calculates the structure and electronic state of a normal molecule or a cluster obtained by cutting out a part of the crystal. Has been used for. However, these methods require a huge amount of calculation time because the Schrodinger equation is solved for each electron in the system. Therefore, it is practically impossible to incorporate these first-principles methods into methods that require hundreds of thousands of calculation steps, such as the molecular dynamics method with empirical potential. On the other hand, in the field of molecular dynamics,
In 1991, AKRappe and WA Goddard "Charge Equilli
In a paper entitled "bration for Molecular Dynamics Simulations", we proposed a new charge calculation method for polyatomic molecules. The present inventors further extend the charge calculation method for polyatomic molecules proposed by ARRappe et al. and, significantly modified to be applicable to various elements and crystal system, Si or SiO ₂ which is generally used as a material for a semiconductor device, of course, be applied strength and dielectric such PbTiO3 and PZT, to the high dielectric material A design simulation method was constructed.First, the method of Rappe et al. Is described below.

First, Rappe et al.
Was considered as a function of ionization potential and electron affinity. That is, in a certain atom A, an electrically neutral state is indicated by A0. Considering up to the second term of Taylor expansion,

[Equation 44] The energy when the atom is electrically neutral is (3031)
It is obtained by substituting 0 for QA in the expression.

E _A (0) = E _A0 (3041) The energy when the atom is +1 charged is (3031)
It is obtained by substituting 1 into QA, and the energy when ionized to -1 valence is obtained by substituting -1.

[0118]

[Equation 45] The difference between (3042) and (3041) is the ionization potential IP
And the difference between (3043) and (3041) is the electron affinity EA.

[0119]

[Equation 46] The reason why the EA in the formula (3045) has a minus sign is that considering that taking electrons is a positive job, adding electrons is a negative job. Taking the difference between (3044) and (3045) gives the electronegativity χA0.

[0120]

[Equation 47] Further, if the sum of (3044) and (3045) is taken, the following equation holds.

[0121]

[Equation 48] Here, consider the physical meaning of the expression (3052). Figure 1 shows a schematic diagram of Coulomb interactions during ionization in the case of hydrogen atoms.
2 to FIG. As you can see, the difference between IP and EA is Coulomb interaction b between electrons. Thus, when the Coulomb repulsion between two electrons in φA trajectory denoted JAA0, it will be referred to as ^{_{IP-EA≡J AA 0 ... (3053}} ) JAA0 and Idempotential. Since the orbital shape changes by adding electrons, JAA calculated by Hartree-Fock's first-principles calculation is slightly different from the value of equation (3053), but here it is approximated by equation (3053). .

Consider the relationship between the Idempotential JAA0 and the atomic diameter RA0. Considering the simple atomic structure model shown in FIG. 15, Coulomb interaction JAA0 between electrons is as follows.

[0123]

[Equation 49] Here, JAA0 is in eV units and R is in Å units. (3053)
In order to confirm the validity of the expression, the relationship between the J value and R is
Summarized in 1. As shown in Table 1, it was found that J is almost inversely proportional to the atom size, and the diameter of the atom obtained from the equation (3054) is approximately equal to the homopolar bond distance.

By using the equations (3051) and (3053), the energy (3031) of the isolated atom can be written by the following equation.

[0125]

[Equation 50] In the case of polyatomic molecules, it is also necessary to consider electrostatic energy between atoms. The electrostatic energy can be expressed by the following equation.

[0126]

(Equation 51) Here, JAB is a Coulomb interaction when there is a one-electron charge at the center of atoms A and B. JAB depends on the distance between AB.
Using equation (3055), the total electrostatic energy is (305
It can be expressed by equation 7).

[0127]

[Equation 52] The atomic scale chemical potential χI (Q1, ..., QN) can be expressed by the following equation.

[0128]

(Equation 53) From the charge balance condition (3059), (N-1) simultaneous equations are obtained.

Χ ₁ = χ ₂ = ... …… = χ _N (3059) That is,

(Equation 54) Becomes Including the condition (3062) for total charge,
Solve N simultaneous equations.

[0130]

[Equation 55] When (3061) and (3062) are displayed in a matrix, they can be expressed by the following formula.

[0131]

[Equation 56] In solving the equation (3063), the electric charge that each ion can take has an upper limit and a lower limit. Therefore, if the electric charge is out of the range, the electric charge is fixed to the boundary value. The equation (3061) is separated into a charge fixed at the boundary value and a charge not fixed, and (3)
063) is modified to solve the matrix and converge. Then, the variation in charge can be calculated.

A problem in this method is JAB's decision. When the distance RAB between atoms A and B is large

[Equation 57] However, it is necessary to consider the repulsive force at the distance where the electron clouds overlap. Here, the electron density is approximated by the term of one Slater orbit. That is, in an atom whose outer valence orbital is ns, np, or nd, in any case, it is assumed to be a standardized ns Slater orbital of the following type.

[0133]

[Equation 58] When the average size of atoms is calculated using the equation (3065),

[Equation 59] Becomes The valence orbital index ζA is selected according to the following relationship.

[0134]

[Equation 60] λ is a parameter for adjustment, and is included to calculate the difference between the average atomic size given by equation (3066) and the covalent bond radius RA of the crystal.

The scaling parameter λ is set by using an alkali metal halide molecule as shown in FIG.
The alkali metal halide molecule is MX, M is an alkali metal (Na, K, Rb, Cs), and X is a halide (Cl, B).
r, I). Representing each charge by QM and QX, (3
062), Q _total = Q _M + Q _X = 0 ∴Q _X = -Q _M (3068) From formula (3061)

[Equation 61] When the expression (3068) is substituted into the expression (3069),

(J _MM −J _MX −J _XM + J _XX ) Q _M = χ ⁰ _x −χ ⁰ _M (3071) By the way, it is clear that J _MX = J _XM (3072)

[Equation 63] The only variable in the equation (3073) is λ. From the calculated Q
The dipole moment μMX is calculated using the equation (3074) and compared with the experimental value to determine λ.

Μ _MX = 4.80324Q _M R _MX (3074) R _MX is an experimental value of the bond distance. In the expression (3074), Q is e unit, R is Å, and μ is Debye unit. From the comparison with the experimental value of dipole moment, 0.49
13 has been obtained. By using this λ, the valence orbital index ζ can be obtained from the values in the table of FIG. 17 from the equation (3067).

Rappe et al. Treat hydrogen differently from other elements.

The electronegativity of hydrogen is inconsistent when compared with the experimental values of Mulliken and Pauling et al. On the Pauling scale, hydrogen is more electropositive than carbon and slightly more electronegative than boron, whereas on the Mulliken scale hydrogen is more electronegative than carbon and nitrogen. This electronegativity problem is because the bound hydrogen orbit cannot spread like free hydrogen ion, so the effective electron affinity EA becomes smaller than the atomic value.

Therefore, using EAH as a variable, hydrogen χ0H and J
Redefine 0HH.

[0140]

[Equation 64] From the comparison with the experiment, by the method of least squares, χ ⁰ _H = 4.528eV J ⁰ _HH (0) = 13.8904eV (3076) From the comparison with the charge calculated from the electrostatic potential of the HF wavefunction, not the experiment , Χ ⁰ _H = 4.7174eV J ⁰ _HH (0) = 13.4725eV (3077) may be used.

In addition, a boundary value regarding charges is set.

Clearly, χ _A ⁰ and J _AA ⁰ are not appropriate outside the range where the valence shell is empty or full.
For example, in the case of Li, C and O, the charge is limited to the range shown in the figure below. Outside this range, E _A (Q) = ∞
Take.

It is judged whether or not the electric charge obtained by solving the matrix (3063) is within the range of the inequalities shown in FIGS. 18 to 20, and if it is out of the range, the electric charge is set as the boundary value. Fix it.

The above is the outline of the method of Rappe et al. R
appe et al. apply the method described above to organic molecules.
The present inventors have adopted only the concept of atomic-scale chemical potential proposed by Rappe et al. And newly developed a charge calculation method for individual atoms in a crystal system. The charge calculation method by the present inventors will be described below.

The second term on the right side of the equation (3057) means the sum of electrostatic energy of all bonds. The present inventors have tried to extend from molecules to crystals and incorporate periodic boundary conditions. Then, when calculating the electrostatic energy, it was discovered that it is necessary to consider not only the hundreds to thousands of atoms used in the calculation, but also the electrostatic energy received from the innumerable virtual atoms outside it. That is, I noticed that the electrostatic energy had a small attenuation and that the influence from distant atoms could not be steamed. Therefore, assuming that there are virtually innumerable cells outside the cells used for the calculation, the electrostatic energy from the virtual cells is also taken into consideration.

In this case, the electrostatic energy JAB received from the atoms in the virtual cell of formula (3057) is considered to be only Coulomb energy because RAB is usually a large value of 10 angstroms or more in consideration of the size of the calculation cell. It is good to change. Therefore, JAB = 1 / RAB is treated uniformly. When the virtual cell is represented by ν, the equation (3051) is changed as follows.

[0147]

[Equation 65] Here, ν = 0 means the calculation cell itself, and ν ≠ 0 means a virtual cell. The Ewald method is applied to the third term of the equation (3078). Using the Ewald method, the third term can be calculated as follows.

[0148]

[Equation 66] G is the reciprocal lattice vector, Ω is the volume of the calculation cell, and κ is the parameter that adjusts the convergence of the Ewald calculation. Furthermore, (30
The first term on the right side of the equation (82) does not include the term of ν = 0 and A = B, and the second term does not include the term of G = 0. From these equations, when the atomic scale chemical potential χI is calculated,

[Equation 67] Becomes From the condition of charge balance, N-1 equations are obtained.

Χ ₁ = χ ₂ = χ ₃ = ... = χ _N (3092) (3091) Substituting the equation,

[Equation 68] It solves N simultaneous equations that can be expressed by equation (3093).

When a matrix of CQ = -D is created from the expression (3093), the following expression is obtained.

[0151]

[Equation 69] The present inventors have discovered that when handling crystals such as SiO ₂ or amorphous materials, the equation (3094) should be solved while considering the periodic boundary conditions. Since the formula (3094) is non-linear with respect to the charge Q, the calculation program has a roof for converging Q.

Furthermore, the present inventors have paid attention to modeling of Coulomb integration in actual calculation. JAB is
It was calculated as the two atom Coulomb integral of the Slater function for ζA and ζB. Coulomb integral JAB. Was calculated using the method of Roothern. Roothern et al. Use the label for the Coulomb integral JAB. The Coulomb integral can be expressed by the equation (3102) using the parameter defined by the equation (3101).

[0153]

[Equation 70] here

(71) (α, α ′, β, β ′) = (1 + τ _a ) α (1-τ _a ) α ′ (1 + τ _b ) β (1-τ _b ) β ′ ... (3103). further,

[Equation 72] Here, the unit system of the above formula is the atomic unit system, and the distance is boh.
Hartree is used for r and energy.

Although Roothern et al. Give an analytical expression of Coulomb integration for 1s, 2s, 2p, the present inventors calculated up to Coulomb integration for atomic orbits after 3s, and use Si commonly used in semiconductor devices. The optimum orbital was selected so that SiO ₂ and SiO ₂ could be handled.

The value of JAB has a great influence on the charge calculation. Therefore, it is important to set the Slater orbital and valence orbital index ζ used when calculating J. Therefore, the present inventors
We investigated the value of J when dealing with Si-O and Si-Si.

Rappe et al. “Assume that the polarization charge is on the ns orbit”, but it is unclear which orbit is optimal to consider for various elements. Since the main quantum number n of Si is 3, we calculated and compared the case where the charge of Si was modeled to be in 3s and the case where it was assumed to be in 2s. At this time, the analytical formula of J at n = 3 was separately prepared by the present inventors according to the method of Roothern. For example, Si
The following equation was used for -Si.

[0157]

[Equation 73] FIG. 21 shows the calculation result of Coulomb interaction J between Si and Si. In this case, it can be seen that there is almost no change in the value of J regardless of whether 2s or 3s is used. However, the Coulomb interaction J between Si and O is greatly different when 2s and 3s are used. In Fig. 22, the charge of O is in 2s and Si is 2s or 3s.
The JAB calculation results are shown for the case of modeling. S
When the charge of i is set only to 3s, the value of J has a local minimum value around r = 1.6Å due to the influence of the 3s orbit spreading very much. Since the bond distance of Si-O is just around this,
It was found that J was estimated to be very small, and as a result, both Si and O charges were extremely underestimated. The same was true for Si-H. Therefore, if we dare to assume the 3s orbit as a model of Si charge, the estimation of J would be completely inaccurate. From this, it was found that the slater orbit was used only for considering the spread of the charge, and the 1s or 2s slater orbit was the optimum model. This is also the case for elements with atomic numbers after Si, and it was found that when the extended orbits after 3s are used as a model, the calculation accuracy of Coulomb force becomes very poor when the interatomic distance is short. .

The present inventors also paid attention to the program for the correction for hydrogen. From equation (3031),

[Equation 74] Atomic-scale chemical potential χHl after hydrogen correction
Ask for. First, assuming that the formula of energy of polyatomic molecule corresponding to formula (3057) is M number of hydrogen among N number of total atoms,

[Equation 75] Here, the symbol Σ ′ in the first term, the second term, and the fifth term on the right side of the equation (3114) means that H is excluded when the respective sums are obtained. The atomic-scale chemical potential χHl is

[Equation 76] The fourth term on the right side of equation (3121) can be solved as follows.

[0159]

[Equation 77] From equation (3122), the sum of the fourth and fifth terms of equation (27) is
It can be seen that it means interatomic Coulomb interaction with all atoms other than Hl. Therefore, the chemical potential is

[Equation 78] Using the formula (3123), the matrix (30
63) will be changed. Now, assuming that χ1 is an atom other than hydrogen, substituting equation (3123) into equation (3059),

[Expression 79] Compared with the matrix (3063), hydrogen Jii0

[Equation 80] Will be changed to.

Next, consider hydrogen correction of the interatomic Coulomb potential JAH. The Slater orbital function for hydrogen is (30
From equation 75),

[Equation 81] If hydrogen is included, after obtaining the charge (3126)
Correct JAH using the formula and recalculate the charge. This is repeated until the charge values converge. At this time, I = 1
When the row of is a hydrogen, the present inventors created an algorithm so that the row is replaced with a row of which the hydrogen is not formed to form a matrix.

Regarding the handling of the charge boundary condition, the present inventors also paid attention to the specific method of changing the matrix.

The electric charge QB fixed to the boundary value by the equation (3061)
And consider it separately into charge QC which is not so.

[0163]

(Equation 82) When the equation (3062) is changed to the equation (3131), the elements of the matrix are

[Equation 83] Change the matrix as in the above equation to keep the charge within the boundary conditions. At this time, when the element with I = 1 is set as the boundary value, the rows of the matrix are exchanged and the row is selected so that the row for the element for which the boundary value is not set becomes I = 1.
Enabled to select (pivot).

The present inventors have used the above-described charge calculation method newly developed for crystal systems and amorphous systems, and adding this to the molecular dynamics method to sequentially calculate the movement of individual atoms to obtain a dipole. The moment was calculated every moment and recorded, and this was Fourier-transformed to calculate the IR spectrum. FIG. 23 summarizes the above procedure.

On the other hand, in the gate oxide film used in the DRAM of 16G or later and the tunnel oxide film of the NAND EEPROM, high reliability is an important issue more than ever before. Especially in the EEPROM, since the oxide film is used under the severe condition that the insulating property is maintained while the tunnel current is passed under the high electric field, it is important to search the limit of the insulating material at the atomic / electron level. Has become.

Under these circumstances, various research institutes proceeded with the atomic level search experiment of the Si / SiO ₂ interface structure and the SiO ₂ network structure, the experiment of the SiO ₂ trap mechanism, and the prediction by theoretical calculation. Has been. However, as in the present invention, there has been no conventional method for associating and utilizing the relationship between the network structure and the IR full spectrum at the atomic level. Slightly IR spectrum Si
Only experimental results relating the asymmetric stretching peak of -O and the stress at the oxidation interface have been reported.

Further, as the interpretation of the experimental results, the discussion has been made only from the vibration analysis of the Si—O—Si molecular structure using the “central force model”, and the actual network is greatly simplified, and the amorphous structure is used. Distribution of Si-O bond length peculiar to quality, S
Conventionally, no discussion has been made considering all the distributions of i-O-Si and O-Si-O bond angles. In the present invention, by considering these various distributions peculiar to amorphous, the correlation between the electrical characteristics and the structure was obtained from a viewpoint different from the central force model.

Here, the points that the present inventors have overcome will be added. The present inventors examined the potential when using the classical molecular dynamics method, and further,
Made improvements. The present inventors first tried to use the TTAM potential developed by Tsuneyuki et al.
We have found the disadvantage that i-O asymmetric stretching vibrations appear on the very low wavenumber side. Therefore, the BSK potential developed by van Best et al. Was used. Their potential showed a peak on the higher wave number side than the TTAM potential, but an expansion / contraction peak appeared at a wave number lower than the actually measured value.

The present inventors have used these existing potentials in order to bring this peak value close to the actual measurement.
The condition of constant charge was removed, and the distribution of charges when the Si—O—Si angle or the like was largely displaced was incorporated into classical molecular dynamics. This charge distribution calculation can be ideally calculated using the molecular orbital method, but in reality, it is impossible at present to handle the structure of several hundreds of atoms at a time to calculate the charge distribution. Even using a super computer is impossible in reality.

Therefore, the inventors of the present invention used the A. R. Rapp
e and W. A. The charge distribution calculation was performed using the method proposed by Goddard III. The cited paper is "Charge Equilibration fo
r Molecular Dynamics Simu
"", J. Phys. Chem., 95, 1
991, p. 3358-3363. The present inventors developed mathematical formulas with reference to their papers, corrected errors in the papers, calculated necessary parameters independently, and applied their method to the SiO ₂ system. The outline of the method of Rappe et al. Will be described below.

Energy Ei by individual charges of atom i
Considering (Q) up to the second order of Taylor expansion, the following equation (1
651).

[0172]

[Equation 84] Here, Xi 0 is an electronegativity and is represented by the following formula (165
It is represented by 2).

X _i ⁰ = 1/2 (IP + EA) (1652) Jii is the Coulomb repulsive force between the electrons on the orbit of atom I and is defined by the following formula (1661).

Jii = IP-EA (1661) Here, IP represents ionization potential and EA represents electric affinity. Total energy E of a molecule composed of N atoms
(Q1, ..., QN) is the sum of the interaction electrostatic energy J _ij Qi Qj between energy and atomic individual atoms.
However, J _ij is the Coulomb interaction of unit charge at the center of atom _ij . Therefore, the following equation (1662) is established.

[0175]

[Equation 85] The atomic-scale chemical potential is calculated by the following formula (166
Given in 3).

[0176]

[Equation 86] This equation (1663) is converted into the condition of charge balance as X1 = ... = X
By substituting for N, the following equation (1671) is obtained.

[0177]

[Equation 87] That is, the following formula (1672) is obtained.

[0178]

[Equation 88] Further, the condition of total charge shown in the following formula (1673) is considered.

[0179]

[Equation 89] J _ij is known and is approximately the following equation (1674).

[0180]

[Equation 90] Here, R is the size of the atom when I = j, and the interatomic distance when I ≠ j. Therefore, when the matrix C and the vector d are introduced, the simultaneous equation shown in the following equation (1681) regarding Q is established. That is,

[Equation 91] When defined as, CQ = -d (1682). In the paper, C1i = Qi and Cd = d, but the present inventors modified and used them as described above. Furthermore, necessary parameters such as the size of atoms such as Si—O—H were calculated independently, and the simultaneous equations were solved to obtain the charge distribution. As a result, the peak wave number became much more consistent with the actual measurement.

FIG. 2 shows an example of "family register" calculated by the molecular dynamics simulator II, that is, an example of enormous data showing bond length (position vector r (t)) and bond angle per unit time.
4 shows. This data is input to the film property calculation unit III, where basic film properties (optical, electrical, crystallographic properties) are calculated.

Regarding the process variation calculation unit I, an example of the two-dimensional general-purpose system is disclosed in Japanese Patent Laid-Open No. 3-16.
4904.

The present inventors succeeded in constructing a three-dimensional process simulator this time. This will be described below with reference to FIG. First, we give an overview of the technology and the structure of the simulator. Although the present application is mainly based on the atomic level, it also proposes a new method of merging with another process variation calculation unit I, and it is also possible to provide a new function which is more effective and accurate, and which has never existed before. Shows. The process variation calculation unit (10) of the other party to be combined is, for example, a three-dimensional process simulator Ic and belongs to a so-called fluid model. In particular, the contents of the three-dimensional process simulator Ic newly constructed by the present inventors will be described below. There are two main points. First, the "three-dimensional process simulator" constructed by the present inventors can be connected to an atomic level simulator. Next, in this "three-dimensional process simulator" part, the chemical models and mathematical formulas newly overcome by the present inventors are used. For the first time, they have become possible to merge with the atomic level. As output from the atomic level simulator, for example, viscoelastic coefficient, stress constant, diffusion constant,
C _ij , etc., which are transferred to the fluid general purpose 3D process simulator. In FIG. 8, (2
3), (33), (43), (53), (63), (7
3) and (83). For example, the diffusion constant of the oxidizing agent or the like is transferred to (23). The method of transfer procedure will be described later.

General-purpose three-dimensional process simulator (10)
Side, according to the input data, the process-induced stress in a wide area of several microns level is calculated in (2) and (3), and the result is once passed through the file (22) and the file (32). By transferring to the molecular dynamics simulator II or returning to the general-purpose simulator (10) again and repeating this, more accurate and more functions can be used. As the calculation procedure, for example, at a certain time, the result of the process-induced stress calculation in a somewhat wide area obtained by the process variation calculation unit is temporarily transferred to the molecular dynamics simulator II, and the correction amount of the physical property constant is calculated here. However, the movement of atoms is again observed, and the corrected physical property constants are sent again to the general-purpose simulator, and the results of the process-induced stress in a wide region are advanced by the general-purpose simulator. By repeating this, if a sufficiently large stress is locally generated, it is possible to see the movement of atoms on the molecular dynamics simulator II side and see whether plastic deformation is possible.

Of course, at the time of impurity redistribution calculation and stress calculation in the above-mentioned general-purpose process simulator part, the point defect behavior calculation part (9) shown in FIG. 8 is called. As point defects, vacancy and interstitial Si in the Si substrate are treated. Without properly incorporating the behavior of these point defects in the diffusion of impurities, the redistribution of impurities cannot be calculated accurately.
In the ion implantation calculation portion (4), of course, a large amount of vacancies and interstitial Si are generated due to implantation damage. These are redistributed during the thermal process, at which time they also interact with impurities and must be reflected in the calculation of the impurity distribution. Also in the silicide / shape / stress calculation part (7), of course, the behavior of the vacancy and interstitial Si in the Si substrate must be fully considered.

By the way, the other (8) in FIG. Oxidation / diffusion
For the ion implantation and the like, each necessary physical property constant is stored in the process variation calculation unit (10) based on many data. For example, regarding the oxidation of Si, PB1 and PB2 of the formulas (1130) and (1131) given later are the same. The diffusion constant of the oxidant in the SiO2 film is represented by D'of the formula (1531) given later. Generally, as described above, there is generally data.

However, data and the like are not always available during LSI development. For example, consider now the prototype development of a memory cell using a ferroelectric film. CD
Consider the case where PbTiO ₃ is formed using V or a sputtering process. Most of the current general-purpose process simulators do not yet know the diffusion coefficients of various impurities in the PbTiO ₃ film, and their stress constants. Also, the segregation constant of impurities at the interface is not well known. In such a case, first enter the step (8), and immediately, as shown in (84),
The molecular dynamics simulator II is once visited, where the potentials of necessary elements and the like are calculated, based on which the molecular dynamics part moves, and the physical property constants are calculated. After this, it is very similar to the part already described,
These physical property constants are stored in the file from (83). Then, the calculation part proceeds in the same manner as the process variation calculation unit.

Here, with reference to FIG. 9, description will be given in accordance with specific processing. Here, the oxidation process is taken up.
That is, a wafer is charged, the wafer is transferred to the substrate processing step, and then the oxidation step is started.

In this oxidation step, for example, the amount of high-purity oxygen gas is set to Pox, and the dispersion thereof is set to σtox. Further, the temperature raising sequence is set to Temp, the variance thereof is set to σtemp, and the valve sequence is set to tox, σtox, and the like. In particular, the valp sequence refers here, for example, when opening the pulp, the details of which are as described in the examples.

For example, in FIG. 8, the oxidation / shape / stress calculation portion of (2) is accessed, and the partial pressure Pox, the temperature Temp and the like are input therein. Then, the oxide growth calculation and the stress calculation can be performed at the same time. At this time, using (21),
Variance / variation calculation is also performed. Here, σpox and σte immediate are input to efficiently calculate the variation from the central value.

From these calculation results, for example, the growth amount of the oxide film, the shape of the interface between the oxide film and Si, and the amount of impurities are calculated. Furthermore, according to the present inventor, stress is also required. Similarly, from the dispersion / fluctuation, the amount of growth of the oxide film, the shape of the interface between the oxide film and Si, and the amount of impurities are calculated for Pox = Pox + σpox and Temp = Temp σtemp. It And
These are once stored in a memory or file (22) and sent to the molecular dynamics part. In the molecular dynamics part, these are input as a signal a, and the position vector r (t) is calculated.

FIG. 10 shows a process diagram similar to that of FIG. 9, but here it relates to the formation of a ferroelectric film.

By the way, a description will be given here of the diagnostic part of the present invention. FIG. 11 shows a part of the binary tree according to this embodiment. Here, a binary tree focusing on the oxidation process is shown. As shown in the figure, first, in
comparing the SiO ₂ film thickness sent from -situ mode two evening, the SiO ₂ film thickness predicted by calculation. And we discuss when the film thickness is within the allowable range. Then, as the next procedure, the 1R vector curp sent from the inlitu monitor is compared with the 1R curve predicted from the calculation. Then, as shown in the figure, first, the spectra are compared, and a case where the spectra are slightly deviated and a case where the spectra are deviated to the high wave number side will be taken up. Then, as the next procedure, it is considered that the residual stress is largely induced in SiO ₂ . So in this case,
The residual stress can be reduced by slowly lowering the temperature. So, modify the annealed part shown in Figure 9,
It is input in (2) of FIG. On the other hand, the spectra are compared, and if they are within the allowable range, the original process sequence is continued.

With respect to FIG. 8, I think that the data flow is almost understood. Then, I would like to mention a little more about the output of FIG. That is, the execution result is temporarily stored in the file (3)
Is stored in Here, a spectrum corresponding to each process obtained in each process and other optical / electrical characteristic results are stored. Although these will be described in detail later, these can be compared with the amount monitored as actual measurement data.

On the other hand, one of the inventors of the present invention has previously proposed a process simulator system using variational calculation of process variation. That is, JP-A-63-174331 and JP-A-3-164904. These proposals can handle processes such as oxidation / diffusion, ion implantation, chemical vapor deposition (CVD), etching and lithography by fluid handling. For example, in the diffusion process using a fluid model, the redistribution of impurities is obtained while distinguishing the total amount of impurities and the amount of electrically activated impurities while incorporating the electrical effect.

The outline of the simulation system shown in JP-A-63-174331 and JP-A-3-164904 will be described with reference to FIGS. 25 to 49. These proposals deal with process variation as follows. For example, taking the processes of oxidation, ion implantation, and lithography as an example, the average temperature (T) and its allowable variation (δ
T), averaging time (t) and its allowable variation (δt), average pressure (P) and its allowable variation (δP), average mask size (L) and its allowable variation (δL), average dose amount (D) And its allowable variation (δD) are input values.

As shown in FIG. 25, the process simulator is provided with a section for calculating dispersion variation in each subroutine of ion implantation, oxidation / diffusion, lithography and the like. As the calculation procedure, a method of calculating the variation while advancing the calculation for these central values is adopted.

The variational calculation will be briefly described. First, the temperature variation in the diffusion process will be described. In the diffusion process, the diffusion coefficient is D, the internal electric field due to the impurity distribution is ε, and the number of electrically activated impurities is N.
The impurity concentration distribution C is obtained by solving the equation based on the law. A procedure of variational calculation will be described with reference to FIG. Here, the left side of the figure is the normal process (Normal pro
cess).

The nonlinear diffusion equation is the following equation (0071)
Can be represented by

[0200]

(Equation 92) Here, k and k + 1 are times, and the impurity concentration Ck + 1 at the time k + 1 is an unknown number. B on the right side is derived from the second term of the following formula (1072).

∂C / ∂t → (C _{k + 1} −C _k ) / t (1072) Since the implicit method is used here, all coefficient matrices appearing in the above equation use D, the average temperature T, and the like. Is expressed.
In the case of a normal process, the predetermined diffusion time tn + 1-th is appropriately divided and solutions are sequentially obtained up to t1, t2, ..., Tn + 1. Here, at each time, t1, t2,
.., tn + 1, the calculation is performed by linearizing these using the Newton method.

Next, the variation process (Process d
(emitted) will be described. Here, the case where the variation δT occurs in the diffusion temperature is shown. Taking the variation of the above equation and knowing the change amount δCK of the concentration at time k, the variation δCk + 1 at time k + 1 is given by the following equation (10
81).

[0203]

[Equation 93] That is, at the time k + 1, the values D, ψ, N, C, etc. at the time of convergence of the normal process are stored and F ′ is used by using them.
Is prepared and the amount of change with respect to the temperature T is prepared, and δCk + 1 is calculated. Here, a part of the matrix elements for creating F'is shown in FIG. The vacancy model is basically used to derive the diffusion coefficient. In this case, the Fermi level changes due to temperature changes, which changes the concentration of charged vacancies and changes the diffusion coefficient. Therefore, first, with respect to the underlying intrinsic carrier concentration ni, the variation is as shown in the following equation (1082).

[0204]

[Expression 94] δn _i = n _i (1.5 / T + 0.605 / kT ² ) δT (1082) Further, changes in the internal electric field of impurities due to temperature can also be expressed using these. Furthermore, the phenomenon of increasing diffusion coefficient (oxidation enhanced) in an oxidizing atmosphere
The variation of (Diffusion Effect) was also considered. In the oxidation process, a linear oxidation rate constant (Line
ar rate constant and parabolic oxidation rate constant (parabolic rate constant).
nt) also changes, so this was also taken into consideration. These are used to create the previous F '.

Here, a typical diffusion process is taken up, and the state of convergence in one diffusion time is shown in FIG. The horizontal axis of FIG. 27 shows the number of iterations, and the vertical axis shows the residual. FIG. 27 shows the convergence status in the normal process. In the normal process, 3 times Newto
It can be seen that it converges at n loop. The variation matrix F ′ shown above is created using D, ψ, N, etc. at the time of this convergence, and the temperature variation is calculated. As shown in FIG. 27, about 50
It can be seen that it converges in a single round and the calculation amount is 1/10 of the total.

Further, FIGS. 28, 29 and 30 show the accuracy of variation. Fig. 28 shows 1000 ℃ as a normal process.
The distribution of impurities in the case of is expressed over a two-dimensional lattice. FIG. 29 shows a variation process, where δT = 1 ° C.
In this case, the calculation is based on FIG. On the other hand, FIG. 30 shows 1001 ° C. as a normal process.

From these figures, it can be seen that there is a relationship of the value of FIG. 30 = the value of FIG. 28 + the value of FIG. Further, when the values in FIG. 30 are plotted again, it becomes as shown in FIG. 31, and it can be seen that when the temperature shifts, the impurity concentration tends to decrease on the high concentration side and increases on the low concentration side. .

There are many processes for producing a semiconductor device, and each process has a plurality of parameters. For example, in the CMOS process
(I) Ion implantation into the N-well region, (ii) P-w
ion implantation into the well region, (iii) patterning of Si ₃ N ₄ for LOCOS process, (iv) LOCOS process,
(V) n + ion implantation, (iv) p + ion implantation, (vii)
It consists of an annealing process. When considering process fluctuations, it is necessary to deal with the fluctuations (or frequencies) of these numerous parameters.

In JP-A-63-174331 and JP-A-3-164904, a combination of the frequency of variation and execution is considered. For the sake of simplicity, let's consider three processes of α, β and γ. In each case, the central value of the process parameter is C (Center), the positively displaced value is U (Upper), and the negatively displaced value is L (L).
and the frequency of those
The definition was made as shown in FIG. 32 from 74331 and JP-A-3-164904.

Considering the steps α and β, the combination is CU, U if the calculation executions such as UU and LL are ignored.
Five combinations of C, CC, LC and CL are possible. When the process of γ is further continued, the seven combinations shown in the figure can be considered. However, in JP-A-63-174331 and JP-A-3-164904, not all of them are treated, but the symmetry of U and L is taken into consideration, and in the case of L, calculation is performed from U.

FIG. 33 shows a part of the search tree of the inference engine, which is the central part of this system. Here, oxidation
We are considering the diffusion process. “Preliminary” in the figure indicates a predicted value by the statistical analysis system, and “actual” indicates an actual measured value.

The allowable value used at this time is a value obtained through all the steps such as oxidation, diffusion and etching. At this time, the ε value is 2
It is obtained by executing the dimensional process statistical analysis system. If the “real” is out of the “pre” ± ε,
As shown in, describe that the oxidation process was abnormal. In this case, the load module of the statistical analysis system is executed depending on how much the temperature or time is deviated. Although the above-mentioned inspection mainly describes the oxidation / diffusion step, the lithography step is also performed. Fluctuations in the lithography process correspond to lateral deviation of the shape of the SiO ₂ film, lateral deviation of the bonding shape, and the like.

The reasoning engine is Common Lis.
Although p is used, this is not only for numerical processing but also for text (L
This is because it is suitable for the processing of ist) and the processing of shapes, and the maintenance of the system is extremely simple. This system not only reads the measured data from the inference engine, but also executes the load module of the statistical analysis system.

According to the simulation systems disclosed in JP-A-63-174331 and JP-A-3-164904, the coefficient of the matrix at the time of convergence of the normal iterative calculation in each step is stored and the linear coefficient is used. By creating a sub-matrix and taking into account the fluctuations of the individual process parameters and the effect of their superposition, about
The variable part can be processed with a calculation amount of 100.
It was also found that, by linking the inference engine based on Common Lisp with the above system, for example, process margins, process abnormality diagnosis, and even optimization can be performed.

In JP-A-63-174331 and JP-A-3-164904, the basic diffusion equation is established for each impurity as shown in the following equation (1111).

[0216]

[Equation 95] Here, x and z are coordinates in the horizontal and vertical (downward) directions of the two-dimensional space, t is time, C is an impurity concentration, and N is an electrically activated impurity amount. Further, D is a diffusion coefficient including concentration dependence and enhanced diffusion effect, and is represented by the following formula (1
112).

[0219]

[Equation 96] In the formula, Di is a basic diffusion coefficient that does not depend on the concentration, and its temperature dependence is expressed by the following formula (1113).

D _i = BB · _exp (−BA / kT) (1113) BB and BA are constants depending on the type of impurities, and k is a Boltzmann constant. FE is given by the following equation (1121), where ni is the intrinsic carrier concentration.

[0219]

(97) And n is _{1/2 [(C D -C A)} + (C D -C A) 2 +4
n _i ² ]. Furthermore, ni is

N _i = 3.87 × 10 ¹⁶ × T ^1.5 × _exp (−0.605 / kT) (1122)

The oxide film growth rate is a parabolic oxidation rate constant R
1 and the linear oxidation rate constant R2 were used to express as shown in the following equations (1123) and (1124), respectively.

R1 = PB3 · _exp (PB1 / kT) ... (1123) R2 = PB4 · _exp (PB2 / kT) ... (1124) When the above parabolic oxidation rate constants are rewritten as B and B / A, At the point of coordinate x, there is a relation shown in the following formula (1125), where Zo is the oxide film thickness, t is the oxidation time, and τ is a constant.

[0222] These are prepared and the variation of each coefficient with respect to the temperature δT is calculated.

Δni = ni. (1.5 / T + 0.605
The variation of n as / kT2) δT is expressed by the following equation (1131)
It becomes as shown in.

[0224]

[Equation 99] Further, the following formula (1132) is established.

[0225]

[Equation 100] The variation of the growth coefficient is expressed by the following equations (1133) and (1134).
Is represented by

[0226]

[Equation 101] By rewriting B = R1 and A = R1 / R2, the following equation (1
141) and (1142).

[0227]

[Equation 102] Then, the variation of t + τ is expressed by the following equation (1143).

[0228]

[Equation 103] The simulation system described above was epoch-making at the time. However, it is no longer applicable to the process control and monitoring technology in recent clustering tools.

Therefore, the present inventors have conducted further studies and completed the present invention. Here, the ferroelectric film of the field effect type MIS device using the perovskite single crystal was optimally designed, and in accordance with this, the device was actually prototyped and evaluated to confirm the completion of the invention.

That is, the eleventh aspect provides that in forming a field effect type MIS element on a semiconductor single crystal substrate, the ferroelectric film of the MIS element is a single crystal, and the semiconductor substrate is In a system including the single crystal ferroelectric film, an evaluation function based on ab initio molecular dynamics theory that expresses the target characteristic is introduced, and the optimum solution of the system under desired main operating conditions is thereby introduced. A method of manufacturing a field effect type MIS semiconductor device, characterized in that the single crystal ferroelectric film and the substrate are designed and arranged.

Also provided by claim 12 are:
When forming a field effect MIS element on a semiconductor single crystal substrate, the ferroelectric film of the MIS element is made of single crystal,
Moreover, in a system including the semiconductor substrate and the single crystal ferroelectric film, ab initio expressing the target characteristic of the system.
Taking the free energy of the whole system as an evaluation function based on the theory of molecular dynamics, the single crystal ferroelectric film is arranged so that the free energy is minimized under the main operating conditions. It is a method of manufacturing a MIS semiconductor device.

Further, according to a thirteenth aspect, in forming a field effect type MIS element on a Si semiconductor single crystal substrate, the gate insulating oxide film of the MIS element is made single crystal, and the semiconductor is In a system including a substrate and the single crystal ferroelectric film, a target characteristic of the system is expressed a
The free energy of the entire system is taken as an evaluation function based on the binitio molecular dynamics theory, and the single crystal ferroelectric film is arranged so that the free energy is minimized under the main operating condition. Field effect SiMI
This is a method for manufacturing an S semiconductor device.

Further, according to a fourteenth aspect, in forming a field effect type MIS element on a Si semiconductor single crystal substrate, the gate insulating oxide film of the MIS element is a single crystal, and the semiconductor is In a system including a substrate and the single crystal ferroelectric film, a target characteristic of the system is expressed a
Taking the free energy of the entire system as an evaluation function based on the binitio molecular dynamics theory, the single crystal ferroelectric film is arranged so that the free energy is minimized under the main operating conditions, and the single crystal A method of manufacturing a field effect type SiMIS semiconductor device, characterized in that the concentration of oxygen defects in a ferroelectric film is suppressed to 0.01% or less.

Further, the present invention provides that in forming a field effect MIS element on a Si semiconductor single crystal substrate, the gate insulating oxide film of the MIS element is made of single crystal, and the semiconductor is In a system including a substrate and the single crystal ferroelectric film, a target characteristic of the system is expressed a
Taking the free energy of the whole system as an evaluation function based on the binitio molecular dynamics theory, and arranging the single crystal ferroelectric film so that the free energy of this system is minimized under the main operating conditions, A method for manufacturing a field effect type SiMIS semiconductor device, characterized in that the local unevenness of energy is set within a deviation value 3σ, and the oxygen defect concentration of the single crystal ferroelectric film is suppressed to 0.01% or less. .

The optimization procedure and trial production results will be described below.

Here, as the single crystal ferroelectric film, for example, PZT-based perovskite is taken as an example, and the Si surface is, for example, the (100) plane. This substrate surface may be, for example, Cu having an FCC structure as a metal electrode. The most important problem is the single crystal ferroelectric film PZT-based perovskite and Si.
It can be replaced by the question of what kind of crystal positional relationship should be placed with the (100) plane. Here, the total energy of the system is taken as the evaluation function which is the essence of the present invention. Moreover, the Si / PZT interface was taken as this system. The present inventors previously created an energy expression formula as a strict evaluation function, and by using this,
Perovskite as PZT under applied conditions, and Si
It has been found that the positional relationship of (100) is optimally set as follows. That is, the midline of the angle between the a1 axis and the a2 axis in the perovskite C4v structure representation is Si (10
0) plane to the [110] axis direction and the c-axis is tilted so that the positions of the first Si and the third Si are S on the Si (100) plane.
It has been found that matching i is the optimal value.
At the same time, it was also shown that oxygen deficiency within 0.01% in PZT was within the allowable range. FIG. 34 shows a conceptual diagram of the present invention thus prepared.

Details will be described with reference to FIG. Figure 34a
The upper half shows a top view of a single lattice on the Si side,
The lower half is a sketch of the PZT plane. Also, FIG. 34c
Point out that it is best to arrange the lines connecting the perovskites in parallel with the [110] Si direction.
At this time, as will be described later in detail, in order for the straight line connecting with silicon to be parallel to [110] Si, the C axis of the perovskite needs to be inclined. By doing so, according to the calculation by the present inventors, it is possible that the total energy is further reduced by 10% as compared with the case where the lines connecting the perovskite and Si are simply arranged in parallel and the C axis is not tilted. I found it. On top of this, PZT
It has been found that if the content is within 0.01%, the original characteristics of the single crystal perovskite are sufficiently exhibited.

The problem of interface energy pointed out here has many influences on device characteristics, as will be described below. Energy due to the mismatch generated at the interface between the ferroelectric film and the substrate facilitates the formation of unpaired atoms,
This promotes the formation of interface states. This interface state hinders the conduction phenomenon and reduces reliability. In addition, defects such as abnormal bands are formed in the film, and there is a possibility of becoming a level. Here, the allowable range of these defects is shown for the first time.

Comparison was also made between the MIS characteristics produced by the present invention and the conventional example. Here, in particular, the CV method is used as a method capable of effectively evaluating the interface state. The characteristics according to the invention also accurately show the variation of the data. According to the present invention, the value of the interface state is lower than that of the conventional example.
It can be seen that it has decreased to about 1/10. Further, the mobility of the MIS element according to the present invention was improved by about 12% when evaluated under the same electric field condition. Also, the reliability has been improved.

The present inventors will explain in more detail the procedure that led to the invention of a new axial relationship that has never existed before. FIG.
Is a completely new strict abinit created by the inventors
It is a structural diagram of io / molecular dynamics simulator. Fig. 36
The upper half is the abinitio part and the lower half is the part of the molecular dynamics simulator. For the first time using this simulator, the inventors have shown what the structure of, for example, a perovskite looks like. FIG. 37 is a perspective view of a perovskite for which the present inventors were able to show a strict structure for the first time. The large balls in the figure indicate oxygen,
The small balls indicate Pb. Here, the Z axis indicates the C axis of the perovskite, and as can be seen from the figure, which surface of the perovskite and which surface of Si are skillfully joined is such a complicated surface that it is extremely inferred. I understand that it is difficult.

As an accuracy check of the above molecular dynamics simulator, attention was paid to the oxide film, and the work of comparing the spectroscopic data was advanced. That is, the rocking mode, stretching mode, etc. of O-Si-O and Si-O-Si were read from the positions of Si atoms and oxygen atoms at every moment, and the natural frequency was calculated from these values. Comparing this with the measured value, 490 cm-1 and 440 cm-1 and 1100 cm-1, respectively.
It is 1111 cm-1, which is a reasonable range. It was also found that the strength of the O-Si-O mode increases when pressure is applied, which is in good agreement with the measured value.

While proceeding with these verification works, the present inventors also calculated the case where a very large compressive stress was applied to PZT. It was also found that the interface having such a structure has a gain of total energy.

The areas used for the calculation will be discussed in some detail. In the calculation, tens of atomic layers are taken in the depth direction on the Si substrate side, tens of atomic layers are taken in the depth direction on the surface, and tens of atomic layers are taken in the left-right direction, and the calculation area is taken. It was Also in the perovskite, the vertical and horizontal heights are several tens of atomic layers each.

It is this method that allows us to find a new positional relationship with a method that has never existed before. The present inventors have developed a strict integration formula of ionic bond potential for a region including a ferroelectric film and a Si single crystal.
The conventional method lacks rigor and cannot accurately determine the energy of the system. The present inventor will sequentially clarify the conventional defects and inaccurate places in this specification.

According to the present inventor, the following can be explained. The present inventors constructed a new simulator which has never existed in the past and, as an example, calculated a single crystal of a ferroelectric film, specifically, perovskite.

The simulator created by the inventor of the present invention will be described below in comparison with the conventional one. For example, first, let us describe the oxide film. Si
Three types of potentials between -O, between O-O, and between Si-Si must be expressed exactly. Between Si-O and O-
Actually, the total potential between O and between Si and Si is extremely difficult, and the cumulative amount of all three types diverges infinitely due to the term of the distance r. In addition, the calculation itself cannot be terminated because the Coulomb potential skirts far, and in the slightly more exact calculation, the conventional method used Ewald's method. However, as will be described later, the calculation lacks rigor.

The method used by the present inventors for actual placement is shown below. Here, the simulator that the present inventors have created for the first time is shown below.

The present inventors have diligently studied and constructed a new unique ab initio molecular dynamics system in accordance with computational physics. This new system made it possible to make new attempts. The details created by the present inventors will be described below.

Although SiO ₂ is used as an example here, it was also confirmed that other materials can be used by simply replacing the element numbers and the like.

FIG. 38 shows the overall configuration of the new simulator proposed by the present inventors. We proposed a new method that synthesizes the molecular dynamics part (MD) and the density functional theory (DFT). First, the use of the molecular dynamics method and the density functional method by the present inventors will be shown below with reference to FIG.

As shown in FIG. 36, the molecular dynamics (MD) part occupies the lower half of the flowchart, and the density functional part (DFT) occupies the upper half. Molecular dynamics method is zero K
It deals with a rather large system at a so-called finite temperature, and uses the potential obtained in advance from the first principle as the potential between atoms (ions). On the other hand, density functional (DFT)
Was used to find the ground state of a small system.
And, naturally, the temperature of the DFT portion was set to 0K. The potential part was obtained from the ground state of the electronic system. Specifically, the wave function (KS equation) was solved and the local density functional was used together.

The calculation procedure of the simulator by the present inventors will be described by taking Si and O as an example. In the DFT part, the total energy of the ground state of the electron system without degeneracy is calculated and set.

[Outside 1] ].

That is,

[Equation 104]

[Outside 2] Interaction energy, the third term is the kinetic energy when there is no electron-electron interaction, assuming the one-body approximation of electrons,

[Equation 105] Is written. The fourth term is the exchange-correlation energy, whose shape depends on the local density approximation. Electric

[Outside 3]

[Equation 106]

[Outside 4] Equation (Kohn-Sham equation)

[Equation 107]

[Outside 5] Must be self-consistant.

[0254]

[Outside 6] Multiply the electron density n by the rugged ε _xc (n) and integrate,

[Equation 108] And Various forms of ε _xc (n) have been proposed such as the formula of Wigner et al. further,

[Outside 7] Instead of using the Coulomb field of as it is, we replace it with a norm-conserving pseudopotential that eliminates the singularity at the nuclear position.
By smoothing in this way, the number of components in the Fourier expansion could be suppressed.

As the total energy, the present inventors use those including Coulomb potential between atoms (ions) for the sake of strictness. That is,

(Equation 109)

[Outside 8] e.g.). From the equation (2071), the first term on the right side is the DFT.
From equation (2065) to the right

[Outside 9] It is. Similar to the case of DFT, that part is calculated from the one-body approximation formula (2065) by Veff

[Outside 10] I put a new idea here. Considering E as potential energy, other virtual kinetic energy

[Equation 110] Is introduced. Therefore, Lagrangian is L = KE (2073) M _I in equation (2073) is the mass of a real atom (ion), but μ and μ are virtual masses, and can be changed according to the purpose as described later. It was Orthonormal condition as a constraint condition

(Equation 111) Therefore, when making the Euler equation, the undetermined multiplier Λ _ik
Is introduced to L

[Equation 112] Take the variation of what is added. as a result

[Equation 113] The temperature T of the kinetic energy in equation (2073) corresponds.
Because of the degree of freedom has a temperature by classical statistical mechanics of energy such as distribution law, it is a real temperature not virtually speaking in relation to the Σ1 / 2M _I R _I ^2. Further, ν _i ν can be suppressed or enhanced depending on the magnitudes of μ and μ.

For example, if μ → M _I is changed from high temperature to T → 0, ψ _i is changed while R _I is stopped and the ground state of the electronic system is obtained. At this time, the left side of equation (2076) becomes 0, and the right side of equation (2071) is

[Equation 114] _Therefore, {Λ _ik } that appropriately forms a linear combination of ψ _i is diagonalized to obtain the KS equation (2064). That is, DF
The same result as solving the KS equation in T is simulated anneal
You got it with aling. However, the dynamical simulation annealing (DS) is performed because the motion of the virtual dynamical system is followed regardless of the Monte Carlo method to create the temperature T state.
Called A).

When integrating the equation (2076), it is not always necessary to treat ψ _i as a variable as it is, but it is possible to treat a expansion coefficient (that is, a Fourier coefficient) obtained by expanding it with a plane wave as a variable. In order to prevent the number of expansion coefficients from unnecessarily increasing, it is necessary to eliminate the singularity at the nucleus position by pseudopotential as in the case of DFT. In particular, when the periodicity of crystals can be used, the linear combination of ψ _i that makes {Λ _ik } diagonal is described as φ _nk .
k is a wave number vector in the Brillouin band, and n is a band index.

Φ _nk is

[Equation 115] Should meet. λ _nk is an eigenvalue of the matrix {Λ _ik } and is an energy level. Expand φ _nk (n, t)

[Equation 116]

[Outside 11] (Not satisfied) also develop the effective potential V _eff

[Expression 117] Substituting equation (2075) and equation (2076) into equation (2074)

[Equation 118]

[Outside 12] Remove it because it will be common

[Equation 119] Get. Since this formula analytically integrates the vibration part, Δt is set larger than that of the purely numerical method.
The amount of calculation can be reduced.

The Coulomb potential is as follows. The present inventor has found that when it is disassembled, it is divided into four sections. In particular, it was confirmed that there was a new fourth term, although there were only three terms in the past. These methods will be sequentially described. As shown below, the present inventors first started to solve the strict equation by using a basic equation in which the dielectric constant is taken into consideration.
First is the basic formula.

[0260]

[Equation 120] The Coulomb potential is different between the conductor ε = ∞ and the vacuum ε = 1, L is one edge of the (cubic) unit crystal, Σ is in the unit crystal, and Zi, oxygen i and the charge and position of the i-th particle. is there. This is because the charge inside the sphere causes polarization on the inner surface of the sphere. There is a layer of dipoles on the inner surface of the sphere that is not a conductor, but the last term in the above equation just works to counteract that effect. Since Ewald's method gives the one on the left side with ε = ∞, the last term of the above equation must be added to obtain the value in the vacuum. Only the results are listed. Even in this part, the conventional derivation is often incorrect.

Coulomb potential

[Equation 121]

[Outside 13] , And set as follows.

[0262]

[Equation 122] It is expressed as Where n _x , n _y , and _nz are x,
y, a vector in the z direction of the ridge, n _x, n _y, (when the bulk crystal) n _z are each an integer ranging + ∞ from -∞. The Σ ′ 'means that j = i is excluded when n = o.

Here, by introducing a new F function,

[Equation 123]

[Outside 14] Found out.

The present inventors further modified this as follows,

[Equation 124] The indices m _x , m _y , and m _z are integers ranging from −∞ to ∞, respectively. When =,

[Outside 15]

[Equation 125]

[Outside 16]

[Equation 126] I used.

Accordingly, the initial Coulomb potential is

[Equation 127]

[Outside 17]

[Equation 128] This does not include force because it does not include. It can be further transformed as follows.

[0266]

[Equation 129] Summarizing the results so far, Coulomb potential is
After all,

[Equation 130]

[Outside 18] / L _x , 0,0) + m _y (0,1 / L _y , 0) + m _z (0,0,1 / L _z ). The good point of the above expression is that the original Σ term attenuates only in the order of reciprocal, whereas in Φ (1), the Σ term is er.
The factor of fc is that the term of Σ in Φ (2) is rapidly attenuated by the factor of exp. Since it works in the opposite direction against the slowdown with Φ (3), select an appropriate κ that balances both. This is the result of calculating the contribution of the Coulomb force in the order of closer distance, and the result when the periphery is a conductor. Another term is added if the environment is vacuum. This is a part that has not been particularly omitted.

The rotation point group created by the present inventors is shown below with reference to FIG. The Si atom is the center of the sphere of radius r (0,0,0)
, And place four oxygen atoms at the vertices of the tetrahedron on the sphere. The position of the C atom can be specified by the Euler angle (θ, φ). To express it in the (x, y, z) coordinates, the vector (0, 0, γ) toward the North Pole should first be rotated about the y axis by θ and then about the z axis by φ.

The rotation matrix is

[Equation 131] Here, when one of the four oxygen atoms is placed at the North Pole (0, *), and φ is indefinite, the positions of the other three oxygen atoms are determined by the undetermined angle ψ, Each ((109 °
28 ', ψ), (109 ° 28', ψ + 120 °), (109 ° 28 ', ψ + 240
°). Therefore, their (x, y, z) coordinates are cos θ = −1 / 3, sin θ = 2√ (2) / 3, and φ are ψ, ψ + 120 °, ψ +, respectively, in the equation (2061). 240 ° is t
Acting on (0,0, r),

(Equation 132) Becomes In order to orient the first oxygen atom pointing to the North Pole in (θ, φ), R (φ, θ) in the equation (2061) can be applied as it is. As a result, the present inventors made the following determination. The (x, y, z) coordinates of four oxygen atoms are

[Equation 133] Becomes The (x, y, z) coordinates of four Si atoms in the unit cell are represented by parameters a, b, and c. When compared with FIG. 40, 2a is a cycle common to the vertical and horizontal directions.

I (0, 0, 0) II (a, b, c) III (ab, a + b, 2c) IV (-b, a, 3c) V (0, 0, 4c) That is, the first When the atom is placed at the origin and the z coordinate rises by a constant value c and is orthographically projected on the (x, y) plane, I, II, III, and IV form a square. The 5th atom is the 1st atom of the following unit cell and has a z coordinate of 4c, which is directly above the 1st atom. The arrangement of four oxygen atoms belonging to each Si atom is such that the relative arrangement to the Si atom rotates by 90 ° around the z axis. The movement of Si atoms as described above is added to it. For the (x, y, z) coordinates of the four oxygen atoms around the first Si atom, the equation (2061) may be used as it is, and the (x, y, z,) of the oxygen atoms around the second, third, and fourth Si atoms may be used. The z) coordinate is the rotation of R (0, 90 °) and the translation of (a, b, c) with respect to the oxygen atom around the first atom, the rotation of R (0, 180 °) and (ab, a + b, 2c) translation R (0,270 °) rotation and (-b, a, 3c) translation can be obtained.

In addition, I II1 and III2, II III1 and IV2, I
II IV1 and V2 pairs are the same oxygen atom, respectively, but these give (of course) the same relations as above. The oxygen atoms I1 and II4 are laterally offset by one cycle 2a.

After conducting the above investigation, the actual process was considered and the defect density allowed was investigated. That is, now, for example, 3200 Si and 32
00 Ti was prepared, and 0 to 100 oxygen defects were formed in this Ti. In this case, the following facts have become apparent from the position of the Si atom and the oxygen atom at every moment. First,
The relationship between the dynamic behavior of oxygen atoms and Ti atoms and the PZT structure factor was investigated. When oxygen defects are present in the single crystal PZT, even Ti and oxygen in the distance of several atoms ahead are obviously disturbed as compared with the case of no defect, and moreover,
The mobility of point defects is extremely high. Comparing these dynamic behaviors with the case of single crystal Si, it was found that the case of PZT had an extremely large influence and range.
Such a vigorous movement has a characteristic that PZT has an ionic bond component and is easily polarized, and therefore, it is considered that a local defect triggers and is linked to a potential change. The inventors of the present invention have made detailed investigations by such a procedure and found that the advantages of the single crystal ferroelectric film can be sufficiently obtained if the content is 0.01% or less.

The atomic level calculation of the ferroelectric material by the present inventor is as follows. The present inventors
We proceeded to understand the structure of ZO and PZT at the atomic level and to calculate the structure and film properties. I will report the calculation results.
The target structure used in the atomic level calculation is (1) P
Complete crystal of TO with oxygen and Ti in highly symmetrical positions, (2) Complete crystal of PTO with oxygen deviated by 0.3 A, (3) Oxygen deficiency in (2) above, c-axis One placed in the ridged portion (Fig. 37), (4) Oxygen deficiency in (2),
One placed along the a-axis (Fig. 43), (5) P
For example, T of TO is replaced with Zr. Especially (5)
In, the artificial substitution of Zr into a layer was also carried out. It was found for the first time that the functions of oxygen and Pb were different in PTO, PZO and the like. From these, the trigger of "spontaneous polarization" of the ferroelectric film was grasped. In addition, the relationship between the atomic arrangement, the strain, the defect structure, etc. of the ferroelectric film and the basic film characteristics was investigated.

First, the atomic level calculation of the most basic PbTiO ₃ film will be shown. When tensile strain is applied or when the Pb position is shifted from the equilibrium point in the trial (Fig. 4
4, FIG. 45) and the case of having an oxygen deficiency, it was investigated whether the most basic charge distribution for obtaining the film characteristics can be obtained. Figure 44 is a picture of a unit cell of PbTiO ₃ .
Further, it is a case where Pb is slightly displaced from the center position.
This is to obtain an isodensity map of the charge distribution at this time. Further, FIG. 45 is also a unit cell, but shows a case where the Ti-oxygen bond portion is distorted. From this result, the following points have become clear. (i) The whole contains a large amount of ionic bond components.
In addition, (ii) it was found that most of the electric charges were attracted to oxygen with strong electronegativity, and then there was a charge around Ti. Further, (iii) the charge distribution is symmetric as it is and spontaneous polarization cannot occur, but it is easily understood that the local electric neutrality is broken in consideration of oxygen “deficiency” and positional fluctuation of Ti. It was also confirmed that the tensile strength lowered the dielectric constant.

Further, for example, in the calculation result when it was observed that there was no thermal disturbance at 25 ° C., the following was found. That is, the ideal stable structure of the PbTiO3 film composed of the point group C4v was obtained. As a result, it was found that even if the oxygen was deviated by a very small amount, it returned to the original point group C4v with stable energy. Is there a more stable equilibrium point, as is usually said, even when the oxygen position is shifted by 0.3Å, which is larger? I investigated this. It is becoming more stable towards a position very slightly away from the starting point of 0.3Å. It is also said that at this position around this point, spontaneous polarization occurs due to this shift. It was found for the first time that there were at least two stable points, one at the position of simple cubic crystal and the other at a position slightly deviated.

The origin of polarization is locally "neutral" in the amount of charge.
Can also be replaced by the fact that the atomic position that breaks down is a stable point in terms of free energy. In the above calculation,
It was found that the presence of oxygen was extremely large. Therefore, it can be seen that oxygen deficiency could possibly be a large trigger both crystallographically and electronically. That is, (1) The charge distribution seems to be localized in each of the three types of atoms.

(2) According to the order in the periodic table, O, Ti and Pb are (1s2,2s2,2p4) and (1s2,2s2,2), respectively.
p6,3s2,3p6,3d2), (1s2,2s2,2p6,3s2,3p6,3d10,4s2,4p6,
4d10,4f14,5s2,5p6,5d10,6s2,6p2). Reflecting this, it can be seen that the charge distribution region of Pb is particularly large.

(3) From a closer look, it can be seen that the charge distribution of Ti is modulated by the immediately adjacent O.

In order to suppress this, F, Cl, Br, I
It turns out that there is a way to successfully replace the halogen elements such as.

Example 2 Although Si.PZT is described in the above example, spinel MgAl is used as a single crystal ferroelectric film.
_{A 2} O ₄ film, cerium oxide CeO ₂ , or strontium titanate SrTiO ₃ can also be used. Further, aluminum oxide Al ₂ O ₃ , zirconium oxide ZrO ₂ , yttrium oxide Y ₂ O ₃ , yttrium-stabilized zirconium Y
SZ, PrO ₂ , calcium fluoride CaF 2, etc., and a laminated film thereof can also be used. Further, the modification in which the silicon PZT is sandwiched at the interface between the single crystal ferroelectric film and the underlying silicon wafer or underlying electrode can be applied to a structure other than this example.

Further, it is important to reduce the free energy under voltage, and if attention is paid only to the distortion energy,
At that time, it does not necessarily have to be the lowest. What we have shown about perovskites is also based on ideas.

An actual measurement system according to the present invention will also be described. FIG. 47 shows a Sawyer-Tower circuit. The circuit of this method will be described. T is a high voltage transformer, Co
Is a known fixed capacitance and CRO is a cathode ray oscilloscope. If Co is chosen to be sufficiently larger than the capacitance of the ferroelectric sample, the secondary voltage of T is almost applied to the sample, but it is divided appropriately by the resistance Rd and added to the horizontal axis of CRO. The abscissa of the horizontal axis of is proportional to the electric field E in the sample. On the other hand, the vertical axis swing of CRO represents the voltage between both terminals of Co, which is equal to DS / Co. That is, the CRO graphic represents the D-E curve, and the absolute values of D and E can be easily obtained from the deflection of the vertical axis. The spontaneous displacement Ds is
Obtained as D when the saturated portion is searched for E = 0, the spontaneous polarization Ps can be Ps = Ds. Observing a lossless paraelectric with this circuit, it becomes a straight line that passes through the origin, but when there is a loss or the ferroelectric has hysteresis. FIG.
In FIG. 8, the DE curve according to the present invention is shown by a solid line, and the result of the conventional PZT is shown by a broken line for comparison. As can be seen from the above, it is also understood that the present invention has a very high dielectric constant.

What is pointed out here is that the energy due to the mismatch generated at the interface formed between the ferroelectric film and the substrate is to be minimized at the time of application. Moreover, considering the practical use under these conditions, the allowable range of the quality is also shown. Reducing this interface energy leads to suppression of formation of levels that hinder the conduction phenomenon and improvement of reliability. This structure, according to the present invention, is a new abiniti
o / Predicted by a molecular dynamics simulator. Since the ab initio / molecular dynamics simulator equipped with a desired evaluation function does not use any empirical model or the like, what is pointed out in this system is that all natural phenomena can be predicted completely. Here, an example focusing on the total energy when the system is applied is shown, but it is suitable for designing the dielectric constant to a desired value, and if deterioration suppression is a desired problem, for example. Design becomes possible.

The method for calculating the two-dimensional physical properties has been described above. An example of the method for calculating the three-dimensional physical properties has already been described with reference to FIG.

Furthermore, let us show how the actual comparison and diagnosis are performed. That is, the simulation result is transferred to the scheduler (13). At this time, on the other hand, monitor data is sent from the actual experimental system every moment. Then, while these are stored in the file (4), the data transfer from the files (4, 3) is synchronized by the scheduler. And
It is transferred to the comparison / verification / diagnosis block in the comparison / diagnosis verification unit IV.

FIG. 49 is a conceptual diagram showing how an in-situ monitor measuring device is added in the clustering tool used in the present invention. FIG.
As shown in FIG. 9, the monitored amount is, for example, an input gas composition analysis, an XPS signal, an FT-IR signal, and the like, which are average values in a predetermined time space and are only a part of information. We provide a way to contrast these limited measurands with the sufficient information gained from the calculations on the other hand. These are displayed as a graphic image on the monitor as shown in FIG.

FIG. 51 shows an ab initio according to the present invention.
It shows a temperature sequence as an example of input data to the molecular dynamics simulation system. FIG. 52 shows the concept of atomic level execution prediction by the abinitio molecular dynamics simulation system according to the present invention. Thus, for example, at the atomic level, the reaction on the substrate and the dissociation of molecules are described.

FIG. 53 is an ab initio according to the present invention.
It is a visualization figure of the atomic level execution prediction obtained by the molecular dynamics simulation system. This is amorphous S
It is the atomic level view of the state of being annealed at a high temperature near 600 ° C. after i is deposited.

The ab initio molecular dynamics simulation system of the present invention can be applied not only to Si but also to an oxide film. That is, FIG.
An execution example of SiO ₂ of the ab initio molecular dynamics simulation system according to the present invention will be shown. Here are sought every moment is SiO-Si angle of SiO _2, from now, Fourier transform, obtaining the natural frequency, by superimposing them, spectrum of FT-IR is determined only by calculation . By using this prediction result, it is sufficient to compare and verify the actual measurement value with the hourly data.

FIG. 55 shows an example of execution at the atomic level using the simulation system of the present invention when the input data shown in FIG. 51 is used. From FIG. 55, it can be seen that it is being executed as expected.

FIGS. 56 to 69 show changes in Si-Si distance, bond angle, time and temperature from Si film formation to annealing and solid phase growth when the input data of FIG. 51 is used. It is the one that has been closed.

In the system shown in FIG. 1, the basic film characteristic calculation is performed by the film characteristic calculation unit III. However, as in the system shown in FIG. 70, the special film characteristic calculation unit III may be omitted and the required characteristics may be calculated inside the process variation calculation unit I. They can be appropriately modified according to design requirements.

As described above, according to the present invention, it is possible to significantly reduce the hardware and the amount of calculation as compared with the conventional one. For example, regarding the calculation amount, in the process simulator TOPAZ, as shown in FIG. 98, the memory is about 24 Mbytes and the execution time is several minutes. In contrast, the famous process simulator SUPREM
In the case of 4, the memory is 200
Mbytes or more, execution time can reach a day with clay. Furthermore, when it comes to molecular dynamics and ab-initio calculation, the memory and the execution time have reached several times.

In such a situation, molecular dynamics and ab-
The cost of a computer will be significantly high for the prediction of the dynamic behavior of atoms and the prediction of various optical characteristics by fully using the initio tool. The present inventors have found a method for overcoming these disadvantages by making full use of new mathematical techniques, and have also found a new application technique of a simulator by this technique.

Conventionally, there are some that add some procedures to the above. That is the flow sheet shown in FIG. 99, and shows the main part of the flow of the density functional in the conventional molecular dynamics simulator. That is, in this method, the virtual mass is given and the Lagrangian of the entire system is used. This corresponds to the equation (1223) previously described by the present inventors. This corresponds to the already known Car-Parrinello method. It is said that this method can be calculated at a seemingly high speed. However, still, for example, a Clay company level super computer requires one day and one night for calculation in one case. Therefore, in order to investigate the fluctuations of various parameters, if the calculation is performed once, it cannot be a real-time diagnostic tool.

[0295]

BEST MODE FOR CARRYING OUT THE INVENTION Various embodiments of the present invention will be described below.

Example 1 This example relates to a direct prediction function of an optical spectrum, which is one of the great features of the present invention. Here again, the method devised by the present inventors will be described in detail.

First, the technique relating to the design of highly reliable oxide film and the technique relating to the evaluation thereof will be described a little from the conventional technical level. Known publications include, for example, “The Physics and Technology of Amorphous SiO 2
₂ '' Plenum Publishing Corporation, ISBN 0-306-42929
-There is 2. This publication describes IR results and Raman results as optical measurement results of various oxide films. This publication is said to be particularly good for discussion at the atomic level.

P. Of this publication. In 63, Paul McMillan discusses how the Raman signal of an amorphous oxide film changes depending on the process. 71 and 72 are the results measured by Paul McMillian. For example, a curve a in FIG. 71 shows a Raman measurement result of a normal oxide film, and a curve b.
Shows the Raman measurement result of the densified oxide film.

Conventionally, it was considered that the spectra of the two were compared, and that the spectrum was changed by desorption of water due to, for example, densification. Further, FIG. 72 shows the spectra before and after the pressure is applied to the single crystal quartz. Especially, there is a change at 610 and 500 cm-1. However, as can be seen from this, at first glance, the amount of change in the spectrum is extremely small, and it requires a great deal of skill as to where to focus and whether to regard it as a significant difference. Furthermore, more complicated experience is required to identify the amount of change in the substance from the amount of change in the spectrum.

In such a situation, in the development of new materials which will be further needed in the future, there are many cases where there is no spectrum identification table, and the interpretation of the spectrum becomes a major obstacle.

Further, p. From 71, Kobayashi et al. Observed changes in Raman spectra of oxide fibers under high stress. A part of the result is shown in FIG. 73. Here, the case where no stress is applied to the SiO ₂ fiber and the case of 3.5
The figure shows the case where a tensile strain of 10% is applied. As can be seen from FIG. 73, when a tensile strain of 3.5% is applied, there is a change in the vicinity of 450 cm ⁻¹ . However, as can be seen here, the amount of change in the spectrum is extremely small, and it requires a great deal of skill as to where to focus and to regard it as a significant difference. Further, further experience is required to identify the amount of change in the substance from the amount of change in the spectrum.

In the above publication, Kobayashi is trying to interpret the spectrum using FIG. As can be seen from FIG. 74, when tensile strain is applied, the SiO ₂ structure 10
It is said that the tetrahedron portion of 9 and 5 degrees is distorted. However, he only brings out a single set of tetrahedral structures, and does not discuss how they are connected to the tetrahedral part next to them. In reality, since it is a certain finite temperature, it must always have a distribution of angles.

As the inventor has repeatedly described in this specification, if the quantitative distribution is not grasped, it is far from the analysis of the spectrum. The same discussion is also developed by Kobayashi in FIG. FIG. 33 (a) shows O
FIG. 6 is a characteristic diagram showing the wavelength of normal mode vibration of a SiO4 tetrahedron as a function of -Si-O bond angle. In addition, FIG.
(B), (c) is a characteristic view showing a decrease in strength due to V4 mode decoupling caused by stress.

FIG. 76 is a p. It is a defect model figure in the oxide film by WBFowler described in 107. FIG. 76 (a) shows a binding model of a part of perfect α-type quartz, where L is a long bond and S is a short bond.
(B) is an E1'-centered model, (c) is an E4 'model,
(D) shows the E2 'model, respectively. (B) ~
The unpaired electron in (d) is represented by e-.

As can be seen from FIG. 76, the above publication discusses the defective structure by preparing two or one 109.5 degree tetrahedral structure. However, as mentioned above,
In practice, it should always have an angular distribution at some finite temperature. Therefore, without some large number of calculations, it is far from spectral analysis.

FIG. 77 is an example of execution by a new molecular dynamics simulator created by the present inventors. Si is 32
And 64 oxygen were used. The reciprocal lattice points were 520 points. The vertical axis represents the rate of change in volume and the horizontal axis represents time. From FIG. 77, it was confirmed that after 1000 ft, it was almost in equilibrium with the external pressure. For example, O-Si-O at this time
The angle, Si-O-Si angle, and O-Si distance are shown in FIG.

As can be seen from FIG. 24, enormous data corresponding to the “family register” of each atom is required. The present inventors have prepared some ideas for the calculation method. For infrared absorption, for example, a vibration mode satisfying the infrared absorption selection rule may be extracted from the time-series data of the vibrational motion of the atom every moment, and the frequency may be calculated. In any case, such an attempt has not been announced yet, and it can be said that this is a new "weapon" for process development.

The present inventors set the time to 0 when the equilibrium is sufficiently reached by using the molecular dynamics system,
After that, the time until time t is set as the sampling time. The time correlation was calculated as follows. That is, It is. Here, M (t) is the polarizability and is represented by the following formula. The integral value of the above equation is a value obtained by integrating from -∞ to ∞.

M (t) = Σe _i · r (t) (1362) Further, <M (t) · M (0)> represents the autocorrelation function of M.

Here, ei is the number of charges at that point, and i is the particle number. R is a vector,
The measurement starting point is not particularly limited, but it is preferable to set it at the corner of the calculation cell.

To obtain the time correlation, s (t1) .s (t2) .s (t3) .s (t4 at each time.
) .S (t5) .s (t6) .s (t7) .s (t8
).

For example, the calculated value of the coordinate position every moment over the entire sampling time, the energy of its motion, and the momentum etc. exist over all the particles. Moreover, the distance between the particles and the angular positional relationship between them are also calculated.

The present inventors have found that some crystalline oxide films,
The calculation was also performed for the amorphous oxide film. A part of the result is shown in FIG. That is, FIG. 78 shows the calculation result of the Si—O distance of single crystal cristobalite as one example.
Further, FIG. 79 shows the calculation result of the Si—O—Si angle when there is a defect. From FIG. 79, it can be seen that the motion is greatly disturbed around the defect. Similar behavior was also obtained by varying the concentrations of these defects.

FIG. 80 shows an example of execution of a new molecular dynamics prepared by the present inventors, which simulates the movement of atoms. In the figure, the large balls indicate oxygen and the small balls indicate Si. 80A shows atoms in the XYZ space, and FIG. 80B shows projections on the XY plane. Although the locus could not be shown here, it was also found that in this single crystal example, the atoms mainly move preferentially in the median direction of the XY axes. This means that a large deviation is shown in the part of M (t) = Σe _i · r (t) (1371) shown in the above equation, and as a result, a spectrum peculiar to a single crystal is shown.

The present inventors performed simulations by changing the measurement temperature, the behavior of individual atoms, the state of natural vibration based on the behavior, and the case where strain was added. 81 and 82 show those examples. That is, FIG. 81 is an X-ray spectrum prediction diagram by the new molecular dynamics simulator according to the present invention, and FIG. 82 is a diagram showing similar stress prediction results.

Example 2 The inventors of the present invention also made it possible to predict the spectrum of an oxide film in the presence of oxygen deficiency and enable comparison with an actual oxide film. Example 2 shows an example in which even the optical properties of the oxygen-deficient oxide film can be calculated.

There is no problem when the ab-initio molecular dynamics method is used, but in the case of the molecular dynamics method using the TTAT potential or the BKS potential, the charge of each atom is treated when oxygen deficiency is treated. The issue is how best to set. this is,
If the strain becomes very large, the bond may be extended or Si-S
This is because i-bonding occurs and a constant charge cannot be assumed. For example, it is expected that when the distribution of charges actually changes due to the elongation of the Si—O bond and the polarization increases, the electrostatic energy of stretching vibration increases and the peak frequency of IR increases.

Therefore, the present inventors have found that Rappe-God.
Dard's equilibrium charge approximation (Charge equilibrati) using the ionization potential and electric affinity that have been tried for alkali and hydrogen.
on approach) was applied to SiO ₂ for the first time, and the electrostatic potential was calculated by incorporating the correction of the charge and the motion and IR of each atom were obtained.

The equilibrium charge approximation is treated assuming that the electrostatic energy of the system is a function of the ionization energy of each atom, the electric affinity of each atom pair, and the charge of each atom. It is a method to solve the simultaneous equations of the number corresponding to the number of atoms under the condition that the potentials are equal, and obtain the charge of each atom.

FIG. 83 shows an IR line spectrum when this method is used for β-quartz and a conventional example. By using this method, the peak frequency was shifted to the low frequency side by about 10%. This value matches the actual measurement with better accuracy.
By using this method, even a-SiO _2, further, it is possible to the structure and IR spectra associated with high accuracy, compared with the actually measured is able to perform effectively.

FIG. 84 is a conceptual diagram showing the method devised by the present inventors in an easy-to-understand manner. The search method is shown in FIG.

The present invention is not limited to the above embodiment. It can be effectively applied to manufacturing other devices such as a liquid phase device and a superconducting device, especially by a clustering manufacturing apparatus. In addition, various modifications can be made without departing from the spirit of the present invention.

Third Embodiment Here, as a third embodiment, details of the direct prediction function of the optical spectrum will be described. Introduction a-SiO
_{A method of making 2} (amorphous SiO ₂ ) on a computer and predicting its optical spectrum will be shown. a-SiO
₂ was prepared by melting single crystal SiO ₂ on a computer and then rapidly cooling it. FIG. 86 shows the temperature setting when forming a-SiO ₂ . The melting may be performed slowly, but here, the temperature was rapidly raised to melt rapidly. After that, the temperature was lowered at the same rate of temperature reduction. FIG. 87 shows changes in pressure at that time, and FIG. 88 shows changes in volume of the calculation region.

In the course of this process, when the empirical molecular dynamics method is used or the ab-initio molecular dynamics method is used, when the temperature is largely changed by temperature increase / decrease, We have found that the handling of external pressure and volume is important. For example, in this embodiment, the external pressure is treated as a function of temperature so that the phase transformation occurs while keeping the volume substantially constant, and the pressure is set to increase as the temperature rises. This is shown in FIGS.

However, if the external pressure is kept constant without paying attention to the volume, a large expansion of the volume occurs as shown in FIG. 89, which is completely different from the actual high-temperature melting SiO ₂ and is remarkably high. It becomes low density SiO ₂ . Furthermore, even if it is in the cooling process or in the solid phase, low density SiO ₂
Is formed. Of course, the potential energy was very high and unstable. If this happens, it will be difficult to compare it with the actual measurement. Therefore, it has been found and executed here that it is desirable to prepare a-SiO ₂ while paying sufficient attention to the external pressure.

FIG. 90 shows fluctuations in potential energy during melting. Since we are paying attention to the external pressure, the potential oscillates while gradually rising. Similarly, the potential energy gradually decreases during quenching. FIG. 91 shows the change over time of two sides of the plane perpendicular to the cell c-axis in the calculation area at this time. The cells can be changed in shape depending on the phase transformation,
-Since the method of Rahman is adopted, it changes due to melting.

Since the quartz is rhombic at first, the angle formed by the two sides is 60 degrees. However, when melted, the crystal structure is broken and internal pressure is evenly applied to each side, so the angle formed by the two sides approaches 90 degrees. It did not return to quartz even when quenched, and seemed to maintain a rectangular parallelepiped at about 90 degrees.

Furthermore, the present inventors have found that there is another point to be noted when forming a-SiO ₂ . It concerns the potentials of Si and O. Ab-initio also has empirical potential (AT
Even in the molecular dynamics method (TA potential or BKS potential), when the SiO ₂ system is treated, the variation of the potential energy due to the variation of the arrangement of the atoms is very large. Ve
This means that it is necessary to obtain the motion of atoms with higher accuracy than the lret algorithm.

FIG. 92 shows the fluctuation of the potential energy and the fluctuation of the kinetic energy, and the sum of the two energies, when the motion of the atom is obtained by the usual Verlet algorithm. If there is no exchange of energy with the outside world,
Total energy should be constant. However, in the usual Verlet's algorithm, the integration of the equation of motion derived from Lagrangian is not performed with sufficient accuracy, and thus energy is not saved in this way. Even if the time interval, which is the integration interval, is set small in order to improve the integration accuracy, the situation remains almost unchanged. Upon examination, this was found to be peculiar to the material due to the extremely steep oxygen potential.

Therefore, in order to handle oxygen accurately, Ver
The let's algorithm was improved, and even higher order acceleration of the motion of individual atoms at different times was obtained, and this was integrated to find the velocity and position. Then, as shown in FIG.
The kinetic energy fluctuates in the opposite phase to the potential energy, and the total energy is kept constant. When such energy is not handled accurately, according to the results of the present inventors, when the movement of atoms is high at a high temperature, liquid SiO ₂ or the like causes poor calculation accuracy, structure is randomly broken, it was not able to create a-SiO _2. Also, the IR spectrum was completely different from the actually measured value.

As described above, the present inventors have found that a-SiO
By Colas their own devising in ₂ of the creation, made it possible to create for the first time a-SiO ₂ on a computer. The a-SiO ₂ thus prepared had a slightly longer Si—O bond length than quartz, and the Si—O length was within a certain range. In molten SiO ₂ , S
The i-O length is S with a bond length smaller than the Si-O length in quartz.
A large amount of i-O was present, but when cooled, this disappeared and only bonds longer than quartz were formed. Examining the coordination number, when the first proximity was 2 A, Si was mostly in the 4-coordination, and 3 and 5 were mixed a little, and the others did not occur. Most of oxygen was two-coordinated.

[0332] Next, over the hydrostatic pressure or uniaxial stress to a-SiO _2, to predict calculate a structural change of a-SiO ₂ in accordance with various stress field, until the potential energy reaches equilibrium, the temperature Was maintained at 900 ° C. to structurally relax SiO ₂ for 2 nanoseconds. The various a-SiO ₂ structures thus formed were individually converted into numerical values such as atomic positions, velocities, accelerations, primary differential values thereof, cell basic vectors, cell basic vector displacement velocities, etc., and stored.

After carrying out the steps up to this point to prepare various a-SiO ₂ structures, IR spectra were obtained for each. The IR spectrum was obtained from the dipole moment. The dipole moment M can be calculated from the following equation M = Σqi · ri (1461). Here, qi is the charge of the atom i, and ri is the position. When obtaining this position, when the system is neutral (Σqi = 0), the value of the dipole moment is the same regardless of where the origin r0 is. that is,

This is because M = .SIGMA.qi (ri-r0) =. SIGMA.qi.ri-.SIGMA.qi.r0 = .SIGMA.qi.ri-r0.SIGMA.qi = .SIGMA.qi.ri ... (1462). If the system is non-neutral, r0
It is understood that the absolute value of the dipole moment changes depending on how the origin is taken, but the Fourier transform only appears in the zero-frequency component and does not affect the others, because the term with is left. The present inventors have found that it is necessary to pay attention to the way of obtaining the atomic position r when obtaining the dipole moment. Since it uses the periodic boundary condition when using the molecular dynamics method, the discontinuous fluctuation of the dipole moment due to the atoms coming out in the −x direction of the cell coming in from the + x direction is considered. To lose.

At the temperature at which IR is actually measured, a-Si
Considering that O ₂ is only vibrating thermally and does not diffuse so much, r used in the calculation of the dipole moment does not come out of the cell in the −x direction without considering the periodic boundary condition. However, the moment was calculated by treating that there was a charge at that position outside the cell. In other words, we decided to handle both atomic positions that consider the periodic boundary conditions and atomic positions that do not consider the periodic boundary conditions when calculating the motion of individual atoms. An example of the time change of the dipole moment thus obtained is shown in FIG. Since it is a three-dimensional calculation, dipole moments in the x, y, and z directions were obtained.

To obtain the IR spectrum, the first step is Fourier transform of this dipole moment to obtain I
A power spectrum that is an R line spectrum is obtained.
One example is shown in FIG. At this stage, the peak position of the IR spectrum can be predicted. That is, about 1000 cm
^A large peak appears at ^-1 and a large peak appears at about 500 cm ^-1 . This result is aS
This coincided with the IR measurement of iO ₂ .

To obtain the IR spectrum, the line spectrum I (ω) is multiplied by the function of ω and the temperature T to obtain the following equation.

[0337]

[Equation 135] α (ω) = (4π 2 / 3hcn) ・ ω ・ (1−exp (−hω / kT)) ・ I (ω) (1463) IR of a-SiO ₂ under various stress fields The stress field, the structure of a-SiO ₂ and the details of the IR spectrum (peak shift, peak shoulder, etc.) were associated with each other. This enabled comparison with actual measurements.

Here, FIG. 96 shows the IR spectrum calculated by applying the molecular dynamics simulator newly developed by the present inventors to α-SiO ₂ , and FIG. 65 shows the newly developed IR spectra by the present inventors. Molecular dynamics simulator aS
Predicted IR when applied to iO ₂ and under external pressure
The spectrum is shown. Further, FIGS. 100 and 101 show the calculation procedure by the present inventors. Note that A and B in the equation of motion of individual atoms in block B in FIG.
Is a term represented by the following formula.

Example 4 By using the system according to the present invention, the magnitude of the induced stress during the process and its components can be predicted, and the detailed relationship between the stress and the physical properties of the material can be grasped. It is possible to predict phenomena such as plastic deformation of the. In this embodiment, the stress problem is taken as an example. We examined how to predict the stress induction and how to reduce it while considering the process fluctuation range.

In the recent fine element process design, it is indispensable to consider the stress when the process progresses. Regarding the stress, for example, the following phenomena occur.

(1) In the buried element isolation method, a groove is formed in a Si substrate and SiO ₂ is buried therein. In this case Si
Due to the difference in the coefficient of thermal expansion between the substrate and SiO ₂ , strain and stress may remain, and crystal defects may be induced depending on the thermal history. Therefore, it is necessary to plan an element structure and manufacturing process that can reduce residual stress as much as possible.

(2) For the purpose of forming a high-quality oxide film, it has become common to use a Si substrate having a low interstitial oxygen concentration. In such a situation, the strength of the Si substrate decreases and it is easy to induce plastic deformation.

(3) In order to maintain the formation of a shallow junction layer and the well region, thermal processes such as annealing and oxidation are performed by R
With the use of TP, the load on the substrate has become even higher.

In such a situation, along with how to suppress the generation of defects, if a stress larger than the specified value or the predicted fluctuation value thereof is observed while monitoring during the process, Immediately clarify the cause,
It is important to show how modifications should be made in future steps. This is even more important with clustering tools.

Furthermore, the type and identification of stress are extremely important. It is important to trace the principal stress and principal surface of maximum stress, but it is also important to focus on the decomposition shear stress. For example, when using a Si (100) substrate, the highest priority sliding surface is the (111) surface as shown in FIG. The [100] direction is the sliding direction on this (111) plane, and it is considered that the most dangerous is when a shearing force acts in this direction.

The simulation system by the present inventors is also capable of outputting these identifications. in this way,
While paying attention to plastic deformation, we will first take up the oxidation of the trench portion, and describe the oxidation shape and stress distribution, as well as the design guidelines for obtaining the prescribed shape. For example, oxidation of the corners of the trench "pointing" or "rounding"
It is said that stress contributes to the phenomenon, but the mechanism is still unknown.

First of all, it is necessary to understand these phenomena in detail and proceed with an accurate process design that suppresses the above-mentioned "sharpness" and promotes "rounding". Further, the shape of the oxide film cannot be accurately grasped unless the induced stress is properly grasped. The problem of stress is not only a problem of plastic deformation, but also an item that is indispensable for grasping the shape and grasping the distribution of impurities under the stress field. In addition, LOCOS of SOI substrates that will be developed in future 1G and later devices
Prediction of the oxide film thickness after the process is particularly susceptible to the stress during the progress of oxidation unlike the ordinary Si substrate, but its quantitative prediction is a more complicated problem.

Against this background, instead of the conventional oxidation-diffusion model, it is necessary to expand the work to grasp the behavior of the oxidizer under stress, understand the oxidation phenomenon under stress, and re-diffuse impurities under stress. I have to. Moreover, it must be calculated at high speed in a short time. The present inventors created a new high-speed general-purpose process simulator part capable of predicting the occurrence of stress and its redistribution, and further calculating the point defect two-dimensional diffusion part under stress, and added this to the atomic level simulator. did. That is, the originally developed high-speed two-dimensional oxidation / diffusion / stress / deformation process simulator was added to the base of the atomic level material design system.

Microscopic Raman and FT-IR were used for actual measurement of in-situ stress observation. The calibration of these measuring instruments was performed in advance. Various physical properties relating to Si and SiO ₂ are obtained by being immediately executed in advance by the atomic level material design system at the time when the process chart is input. This is sent to the high speed two-dimensional oxidation / diffusion / stress / deformation process simulator section.

The inventors first described the phenomenon analysis of "pointing" and "rounding" oxidation which was found by using this new high-speed two-dimensional oxidation / diffusion / stress / deformation process simulator part, and further applied. As an example, the results of a basic examination of the buried element isolation process are shown. When Si is oxidized to SiO ₂ , volume expansion occurs. At this time, SiO ₂ has viscoelasticity and some of the stress is relaxed, but the stress still remains in both Si and SiO ₂ .

FIG. 103 shows such a situation by taking an example of oxidation of the trench corner and examining the situation. First, to reach oxidant concentration C in SiO _2, but going from entering a SiO ₂ film, when the diffusion of the oxidant D'is SiO ₂
In addition to being affected by the stress σ in the film, Si / SiO ₂
At the interface, even when the oxidant reacts with Si, it is affected by residual accumulated stress in the substrate.

As a concrete effect, for example, SiO ₂
The presence of compressive stress σ narrows the atomic spacing of the network and suppresses the diffusion of the oxidant. Also,
If the compressive stress σ accumulates on the surface of the Si substrate,
The Si atom interval is narrowed, the invasion of the oxidant is suppressed, and as a result, the oxidation rate is suppressed.

Now, let the component of the principal stress S in the x direction be σxx,
Further, when the y-direction component is σyy, the relational expressions of the following equations (1531) to (1534) are established regarding the stress on the two-dimensional surface. Where μ is the viscosity coefficient, ν is the Poisson's ratio, and u is the displacement.

[0354]

[Equation 136] By the way, generally, the temperature dependence of the viscosity μ is expressed by the above equation (15
34). Here, VISC. O, VIS
C. E indicates a proportional constant and activation energy, respectively. Therefore, when the temperature dependence of μ is added to the above angular stress equation, it is understood that the components of each principal stress decrease with increasing temperature. Further, the reaction constant κ s ′ when the oxidizing agent reaching the Si / SiO ₂ interface reacts on the Si surface is
Similarly, it can be expressed as the above equation (56). Where κs
Is the value when no stress is present. Where σn
Is σn =-(σxxnx 2 + σyyny 2 + 2σxynx n
y) and στ = − (σxxny 2 + σyynx 2 −2σ
xynx ny). VR is a parameter. From this equation, it can be seen that as the compressive stress increases, the oxidation reaction constant decreases. This is the effect of residual accumulated stress in the substrate when the oxidizer reacts with Si at the Si / SiO ₂ interface.

Furthermore, the diffusion coefficient D'of the oxidant in the oxide film
Can be expressed by the above equation (1531). Here, p is expressed by the following equation. That is, p = -1 / 2 (σxx + σyy). D is a diffusion coefficient in the absence of stress.
The meaning of this equation can be interpreted as follows. That is, the stresses σxx and σyy are such that the tension is positive and the compression is negative. Therefore, it is shown that if there is compressive stress in the film and the substance is in a dense state, the diffusion coefficient becomes small. The above equation expresses this.

FIG. 104 shows the calculation procedure for only this portion found by the present inventors. First, the concentration of the oxidant at each node at the Si / SiO ₂ interface is determined. Then, based on this concentration, a temporary oxidation rate at each node for the first time is calculated. Then, the growth amount u of the SiO ₂ film at each node is obtained. Using this u, the stress S is obtained from the above equation (57). Using this stress, equations (56) and (153)
1) and then recalculate the oxidation rate using equation (1534). Then, again, the growth amount is obtained and the stress is calculated. Then, when the stress is sufficiently converged, a new oxidation time starts to advance again. That is, the Newton method is used.

As can be seen from FIG. 104, the stress calculation during oxidation requires a very long calculation time. Oxidation at 900 ° C. was attempted with a typical trench structure and an initial node of 1300 points. The oxidizing atmosphere is 100% dry oxygen,
The oxidation time was 300 minutes. The calculation time was examined using EWS Sun4. In the case of the conventional technique, it takes about 15900 minutes. This reaches almost 10 days, and as it is, it is extremely unsuitable for practical use.

Regarding the calculation of the Newton's method part, the present inventors have introduced the deviation principle, saved the previous convergence state, and proceeded with the convergence calculation by using the first derivative of this.
With this method, it was possible to reduce the temperature to one thousandth or less of that of the conventional method, taking 900 to 1100 ° C. as an example.
As a result, this simulator has reached the point where it can be used as a real-time tool by optimizing the material property values and increasing the speed of the calculation unit. The inventors further examined the accuracy of the time-dependent value of the oxide film thickness, the shape reproducibility of the trench corners, and the like as the basic accuracy of the simulator.

FIGS. 105 to 77 are graphs showing the results of examining the temperature and time dependence of the grown film thickness using the above simulator. I didn't show it here,
The present inventors also investigated the decompression behavior when the oxygen partial pressure was 10% or even 1%. As a result, when the oxidation rate is slow, the stress accompanying the expansion of the oxide film proceeds while being sufficiently relaxed, so that the stress does not accumulate on the substrate, and therefore a uniform oxide film is formed also on the convex portion. As a result of careful examination by the present inventors, it is effective that the oxygen partial pressure is lowered at the initial stage of oxidation to slow down the oxidation rate and gradually increase the oxygen partial pressure while sufficiently relaxing it. all right.

The present inventors have used this system to investigate why "sharpness" and "rounding" of trench angles occur.
As already mentioned, several factors need to be considered in trench corner oxidation. The following describes five corresponding important elements.

(1) Oxide film viscoelastic model (2) Silicon elastic model (3) Stress dependence model of oxidant diffusion (4) Plane orientation dependence model of oxidation rate (5) Si surface of oxidizer / Si reaction Stress Dependent Model The present inventors surprisingly find that the most dominant of these five models is the fifth model, namely Si / SiO _2.
It was found that it was the stress of the Si substrate at the interface. Here, dry oxidation at 900 ° C. will be taken up.
When the shape of the oxide film and the shape of the Si / SiO ₂ interface at this time are examined carefully, almost no "sharp" shape is observed. 9
After 300 minutes of 00 ° C dry oxidation, 900 ° C
The phenomenon that the trench angle is sharp at the oxidation temperature in the vicinity was analyzed as follows.

That is, no stress is generated so much for up to about 100 minutes. However, due to the SiO ₂ volume expansion along with the oxidation, the compressive stress concentrates in the corner SiO ₂ film and the Si tip side. Especially, in the vicinity of 900 ° C., the viscous flow effect is small and the shape relaxation effect of the oxide film is small. The compressive stress on the Si surface suppresses the oxidation reaction, resulting in sharpening. However, at 1100 ° C., the viscous flow effect is large and the shape relaxation effect of the oxide film is large. Then, there is no accumulation of compressive stress on the Si surface, resulting in rounding. This situation is shown in FIG. As a result of Raman measurement and IR measurement, the predicted values of these stresses could be measured and confirmed with an accuracy within ± 10%.

With the above knowledge prepared, the inventors of the present invention further implemented the above-mentioned simulation system with the STI (Shallow TW) for the purpose of reducing the area of the element isolation region.
It was applied to a part of the basic study of the French Isolation. Here, especially, the basic study of the rounding process of the upper corner of the trench was advanced. The oxidation temperature was 1100 ° C., the oxygen partial pressure was 3%, and it was observed after 100 minutes. As a result, the following can be seen.

(1) Initially, the oxidation of poly Si and the gap in the recessed portion of the buffer oxide film are oxidized. At this time, the corners are slightly sharp.

(2) When the oxidation continues and the gap is oxidized and filled up, the supply of the oxidant is rapidly reduced,
Oxidation from the side becomes rate-determining, and the corners eventually become rounded. The main stress is large near the gap.

On the other hand, 1000 ° C. oxygen partial pressure 50%
The stress distribution during oxidation was investigated. The time is after 30 minutes. As a result, it can be seen that the stress is slightly higher than the stress distribution of 1100 ° C. oxidation.
From the above consideration, to suppress the sharpness of the trench corner,
It is understood that the first condition is that residual stress is not accumulated on the Si surface. The basic measures for this are:
It is conceivable that the oxidation rate is sufficiently suppressed even at a low temperature and the deformation relaxation of the oxide film is sufficiently caused.

The present inventors have further promoted integration with a plastic deformation simulator. As a result, not only oxidation but also T
It was found that the deformation after deposition of the EOS film can also be investigated.

Based on the above results, the present inventors pay particular attention to the decomposition shear stress on the (111) plane, which is a slip plane, and in the two-dimensional simulation system, the value projected on the (110) plane. I saw. In addition, the heat history is 300 degrees /
However, it was found that thermal non-uniformity occurs in the wafer, and if the temperature reaches 600 degrees / minute, dislocations are generated due to stress depending on the part.

The above preliminary knowledge was prepared separately, and the present invention was applied to the manufacturing process of the bipolar CMOS in real time. That is, the comparison between the sequential prediction by simulation and the minimum required in-situ measurement amount in each process, and further, the variation amount and the allowable amount of the factor in each process were examined. Here, in particular, an example in which the variation of the sacrificial oxidation time and the allowable amount of the shape and the like are performed will be described below.

The basic steps are as shown in FIG.
Furthermore, RIE etching for element isolation → oxidation →
There are steps of polycrystal Si embedding → CMP → sacrificial oxidation of the polycrystal Si portion →. There are variations in the factors of each process. In this example, the amount of oxidation was particularly focused on the sacrificial oxidation portion of the polycrystalline Si portion. As the fluctuation range of the oxidation amount, the results of real-time calculation are shown for two cases of zero oxidation amount and a predetermined amount. The amount of oxidation, the oxidized shape on the way thereof, the stress distribution, and further the change of the oxide film shape and the component characteristic of the stress in the subsequent LOCOS process are shown.

FIG. 110 is a distribution chart of σxx at the end of the LOCOS process at 1000 ° C., and FIG.
Distribution diagram of σyy at the end of ℃ and LOCOS process
2 is σxy at the end of the LOCOS process at 1000 ° C.
FIG. 113 is a distribution chart of σxx when the temperature is lowered to 25 ° C. after the LOCOS process is finished, and FIG. 114 is the same as L
After the OCOS process is finished, the distribution diagram of σyy when the temperature is lowered to 25 ° C., and FIG. 115 is the distribution diagram of σxy when the temperature is reduced to 25 ° C. after the completion of the LOCOS process, respectively.

As can be seen from these figures, if there is no sacrificial oxidation on the polycrystalline Si film after CMP, the polycrystalline S
The surface of the i film and the surface of the nitride film are aligned. If the LOCOS process is performed in this state, the oxide film expanded by oxidation pushes from the side of the nitride film. On the other hand, if the sacrificial oxidation is performed to some extent, the nitride film bends well, and the stress relaxes. The result is shown in FIG.

Further, the inventors of the present invention used the above-mentioned simulation system combining the high-speed molecular dynamics simulator using the empirical interatomic potential and the general-purpose process simulator proposed by the present inventors. It was found that the stress remaining in the Si substrate and the oxide film of the MOS element after the process, the interlayer insulating film forming process, and the wiring process has a great influence on the electrical characteristics of the oxide film.

The inventors of the present invention predict the successive stress generation process and the generation site in real time by using the above simulator, and show, for example, what kind of retention characteristics the tunnel oxide film of the EEPROM element part finally has. I also output it. Furthermore, by using this system,
With respect to the allowable width of the retention characteristics of the finally obtained element,
It is also possible to predict what characteristics the device in progress gives. This is the power of the system according to the present invention. Furthermore, during the progress of the process, the system can be successfully used to modify or avoid subsequent processes so that optimum characteristics can be obtained. This example is shown below.

Example 5 This example shows an example in which the present invention is applied to the problem of stress in a silicon thermal oxide film. The stress problem is important for suppressing the generation of crystal defects and dislocations in the silicon substrate,
Not only that, it is also important as a factor affecting the film quality in various materials.

In the silicon thermal oxide film, oxidative stress caused by volume expansion when silicon is oxidized and thermal stress caused by temperature rising / falling occur, which affect the quality of the oxide film. That is, at the atomic level, stress affects the network structure of the silicon oxide film, and S
It may affect the tetrahedron structure, which is the unit structure of the iO ₂ system. It also affects the electrical insulation. Therefore, the present inventors proceeded with the oxidation process while controlling the quality of the oxide film while preparing the knowledge about the IR spectrum, the structure at the atomic level, and the electrical characteristics of the silicon oxide film.

First, the present inventors examined the IR spectrum of the silicon oxide film and the structure and electrical characteristics at the atomic level as a pre-stage. As shown in FIG. 117, a silicon (100) substrate was prepared and an element isolation process was performed to form an element region separated by the element isolation region. Next, the silicon substrate was introduced into an oxidation device to form a gate thermal oxide film as shown in FIG. 118. At this time, the oxidation temperature was variously changed from 800 ° C to 1100 ° C.

After forming a silicon oxide film having a predetermined thickness, the substrate temperature was lowered. At this time, the substrate temperature was lowered in two steps as shown in FIG. The first step is a step of lowering the substrate temperature to 900 ° C., and the substrate temperature was rapidly cooled. The second stage is a step of lowering the temperature from 900 ° C. to room temperature, and no particular attention was paid to the temperature lowering rate.
When the oxidation temperature was set to 900 ° C. or lower, the temperature was lowered only in the second stage.

The setting of the temperature lowering condition was separately determined from the analysis of the viscosity relaxation behavior of the oxide film. It is generally said that the oxide film becomes soft at high temperature and the stress is easily relaxed by viscous flow. Therefore, the present inventors have investigated the stress relaxation behavior due to viscosity in detail. An example of investigating a change in stress relaxation characteristics depending on the temperature decrease rate is shown in FIGS.
1 is shown. In this example, the oxidation temperature was 1000 ° C.

FIG. 120 shows that after performing a typical oxidation process,
When held at 1000 ℃, 600 ℃ / min
It is the result of examining the maximum value of the principal stress in the direction along the surface in the case of rapid cooling at 6 ° C./min. From the results shown in FIG. 120, it was found that the stress was relaxed in about 1 minute in the case of holding at 1000 ° C. and in the case of slow cooling, the viscosity relaxation was slight in the rapid cooling, and only the thermal stress was superposed.

FIG. 121 is a diagram in which stress behavior in the same temperature decreasing process is plotted with the horizontal axis representing temperature. FIG. 121
From this, it is found that stress relaxation is remarkable in the temperature range of 950 ° C. or higher, and the behavior of a viscous body or a viscoelastic body is exhibited, and that stress relaxation hardly occurs at 950 ° C. or lower, and the oxide film exhibits an elastic body behavior. It was

As inferred from such a calculation example,
By using the two-step temperature reduction performed in this example, it is possible to suppress the viscosity relaxation by passing through the temperature region where the viscosity relaxation occurs in the shortest possible time, and by changing the oxidation temperature variously, It is possible to set various thermal stresses caused by the difference in thermal expansion coefficient between the oxide film and the silicon substrate.

Using the sample thus prepared, the tensile stress in the substrate surface direction immediately below the silicon oxide film was measured from the cross section using microscopic Raman spectroscopy. As a result, various thermal stresses were induced in the silicon substrate. Was confirmed. The present inventors also compared the stress-induced measured values in these substrates with the calculation results illustrated in FIG. 122 to infer the stress values induced in the oxide film. In addition, in performing these stress calculations, the present inventors used the material and process prediction simulation system constructed by the present inventors from the molecular level to the general-purpose level. as a result,
The relationship shown in FIG. 123 was obtained between the oxidation temperature and the compressive stress in the oxide film. From FIG. 123, it can be seen that the maximum value of the compressive stress in the oxide film linearly increases as the oxidation temperature rises.

As the next preparation, the IR spectrum was measured for the oxide film having the compressive stress formed as described above. On the other hand, the computer simulation of the oxide film structure and IR spectrum when this compressive stress was applied was performed using the classical molecular dynamics method. FIG. 124 shows a change in the calculated IR vector depending on the presence / absence of the compressive stress.

FIG. 125 shows the result of measuring the wave number shift due to the stress of the Si—O—Si asymmetric stretching peak from the measured value of the IR spectrum. Comparing the measured value and the calculated value, it is found from the calculated IR spectrum that the appearance wave number of the asymmetric expansion / contraction peak of Si—O—Si due to the compressive stress shifts to the lower wave number side than that in the case where no stress is applied. , It turned out that the experiment and the calculation agree well. Therefore, the amorphous structure and IR full spectrum when various stresses were applied were calculated and prepared.

Further, as the next preparation, the relationship between stress and electrical characteristics was examined. Qbd was measured by using the above-mentioned oxide film that produces various compressive stresses as a gate insulating film, and depositing a gate electrode thereon. In order to prevent the viscosity relaxation of the oxide film, all steps such as the deposition of the gate electrode were performed at a temperature of 900 ° C. or lower. Then, the correlation between the oxide film structure predicted from the above-mentioned molecular dynamics calculation and Qbd was examined.

FIG. 126 shows the S of the oxide film under various stress fields.
After calculating the radial distribution function of the i-O atom, S of the first proximity is calculated.
It is the result of having calculated | required the correlation with the frequency of appearance of the bond in which i-O atom spacing is 1.8 angstroms or more, and Qbd. From FIG. 126, when the frequency of appearance in which the Si—O atomic spacing is 1.8 A or more exceeds 5%, Qbd is remarkably reduced, and the Si-0 atomic spacing and the oxide film dielectric breakdown are associated with each other. I found out. The IR spectrum shift amount was about 15 cm ^-1 .

The change in the amorphous network structure when a compressive stress is applied to the oxide film is (1) the expectation that the Si-0 bond length will decrease because the density of the oxide film increases, and (2) From the IR measurement of the thin thermal oxide film and the central model, it is expected that the Si-O-Si bond angle will decrease.
The street has been predicted. However, as shown in FIG. 127, the inventors of the present invention have shown that, as shown in FIG. 127, the increase in the density of the oxide film does not lead to the simple decrease in the Si—O bond length, but rather the spread of the bond length spreads. We have newly found that the proportion of stretched bonds increases. In addition, the correlation with Qbd was also found.

Regarding the physical reason why Qbd and the extended Si—O bond are correlated, the present inventors have not yet clarified the details, but one of the reasons is as follows from a recently reported paper. Guess. The cited paper is "Local deformation in oxide film and hole trap by impurities" reported by Mr. Kaneda in "Research Report on Thin Film and Surface Physics Subcommittee, Formation, Evaluation and Reliability of Ultrathin Silicon Oxide Film" p ． 101-104
It is. Mr. Kaneda is the molecular orbital method program GAUSSIA
According to the calculation using N, "stretched Si-
The presence of O-Si induces the trapping of holes into oxygen at the center. Since the trapped holes are strongly localized in oxygen,
It acts as a trap center for electrons. As a result, most of oxygen captures electrons and returns to its original state. However, a part of them obtains energy of several eV released by recombination of electrons and holes, breaks Si—O bonds that are extended and weakened, and jumps out to the interstitial space, and secondary electrons such as oxygen deficiency are generated. Can generate a specific defect. Is reported.

According to the present inventors, as in this report,
It is expected that stress bonding will increase the number of extended bonds and form oxygen deficiency defects. Further, although various structures can be considered for oxygen deficiency defects, in a structure in which Si atoms at both ends of the oxygen deficiency are largely separated, the electronic state also largely changes, a level is formed in the band gap, and electric conduction It is presumed that the electron distribution will change. It is presumed that the insulating property of the oxide film deteriorates.

By making the above preparations, the IR spectrum of the oxide film in the stress field can be predicted.
Since the correlation between the structure and the electrical property at the atomic level was obtained, the present inventors sought to “correct the manufacturing process”.

The Si wafer which has advanced to the predetermined process is introduced into an oxide film forming apparatus, and is dried at 1050 ° C. containing 5% hydrochloric acid.
An oxide film having a desired film thickness was formed in an O ₂ atmosphere. After forming the oxide film, the wafer was rapidly cooled to room temperature, and the IR spectrum of the formed oxide film was measured. The Si—O asymmetric stretching vibration peak wave number of the IR spectrum was read from the measured value, and the shift amount from the same wave number when stress was not generated was calculated. As a result, a shift of 17 cm ⁻¹ was measured.

From this shift amount, about 0.2 GP is applied to the oxide film.
It was predicted that the compressive stress of a was induced in the oxide film. This shift amount was compared with a shift critical value of 15 cm ^-1 . The thermal process was stopped because it was detected to be above the critical value. Then, as a modification of the manufacturing process, stress relaxation was performed by high temperature annealing.

Further, in order to compare the effects of this example, some wafers were allowed to pass through the process without adding the high temperature annealing process. The annealing time, the annealing temperature, the rate of temperature decrease, etc., which are the parameters of high temperature annealing for stress correction,
Analysis of relaxation time by viscous relaxation prepared in advance, and
It was calculated in consideration of the re-diffusion behavior of the impurity diffusion layer already formed in the substrate. In this case, annealing was performed at 1000 ° C. for 30 seconds.

After the annealing was completed, the IR spectrum was measured again and the shift amount was calculated. Since it was below the critical value, the process was continued as it was. Qbd measurements were made with and without modification of the process. As a result, it was 45 C / cm ² when it was corrected and 5 C / cm ² when it was not corrected, and the effect of the process correction according to the present invention was clearly obtained.

Also, in a wafer having another structure,
A gate oxide film having a desired film thickness was formed in a dry O ₂ atmosphere at 950 ° C. After forming the oxide film, the wafer was rapidly cooled to room temperature, and the IR spectrum of the formed oxide film was measured. As a result, the peak shift amount was 16 cm ⁻¹ . Since this exceeds the critical value, it was necessary to add annealing for stress relaxation, but considering the re-diffusion of the impurity diffusion layer, there was no suitable annealing condition solution. Therefore, it was judged that the process could not be continued for this wafer, and the wafer was discarded.

[0397]

As described above, according to the present invention,
A device that not only solves the problem of the monitor function in clustering but also manufactures a semiconductor device such as MOSLSI,
Especially in single-wafer or clustered manufacturing equipment,
The maximum amount of “in-situ” measurement with limited amount of qualitative content or limited amount of “in-situ” measurement in each device can be expanded to the desired process without test pieces, or It will proceed while making corrections. Further, according to the present invention, (i) the fluctuation of each process element and the superposition phenomenon of the deviation of each element are diagnosed, (ii) the execution of a desired sequence is checked, and (iii) the discovery of a variation factor. Is also possible. Furthermore, the features of this system are (iv) the process variation can be handled mathematically and statistically and efficiently, (v) the addition of a new inference engine, and (vi) the new inverse. The point is that the direction system is added.

[Brief description of drawings]

FIG. 1 is a diagram showing an entire system according to the present invention.

FIG. 2 is a diagram showing that the inclination of the magnitude of the potential with respect to each direction becomes a force in each direction.

FIG. 3 is a diagram showing a new abinitio molecular dynamics simulation system used in the present invention.

FIG. 4 is a diagram showing the effect of entering the action of the first atom (i) and the second atom (j) from the third atom (k) through the second (j).

FIG. 5 is a diagram showing a master-slave relationship of these quantities, where r is a distance between particles and θ is a particle angle.

FIG. 6 is a diagram showing the positioning of a conventional calculation method and the higher-order cyclic partial differential method of the present invention.

FIG. 7 is a diagram illustrating cooperation between a molecular dynamics simulator II and a process variation calculation unit I.

FIG. 8 is a diagram illustrating cooperation between a molecular dynamics simulator II and a process variation calculation unit I.

FIG. 9 is a diagram showing an example of an input sequence of a general-purpose 2 / 3-dimensional process simulator.

FIG. 10 is a diagram showing an example of an input sequence of a general-purpose 2 / 3-dimensional process simulator.

FIG. 11 is a diagram showing a part of a binary tree according to an embodiment of the present invention.

FIG. 12 shows a schematic diagram of Coulomb interaction during ionization.

FIG. 13 shows a schematic diagram of Coulomb interaction during ionization.

FIG. 14 shows a schematic diagram of Coulomb interaction during ionization.

FIG. 15 shows a simple atomic structure model.

FIG. 16 shows an alkali metal halide molecule.

FIG. 17 shows a table of valence orbital indices ζ.

FIG. 18 shows charges obtained by solving a matrix (3063).

FIG. 19 shows charges obtained by solving a matrix (3063).

FIG. 20 shows charges obtained by solving a matrix (3063).

FIG. 21 shows calculation results of Coulomb interaction J between Si and Si.

FIG. 22 shows JAB calculation results in the case where the state where the electric charge of O is 2s and Si is 2s or 3s is modeled.

FIG. 23 shows a summary of procedures for the above-described charge calculation method newly developed for crystalline systems and amorphous systems.

FIG. 24 shows O-Si-O angles, Si-O-Si angles, and O that are execution results of the molecular dynamics simulator according to the present invention.
-The data figure which shows Si distance.

FIG. 25 is a diagram showing each subroutine including variational calculation in a conventional simulation system.

FIG. 26 is a diagram showing a procedure of variational calculation in a conventional simulation system.

FIG. 27 is a characteristic diagram showing convergence states of a conventional process and a variation process.

FIG. 28 is a data diagram showing the accuracy of variation.

FIG. 29 is a data diagram showing the accuracy of variation.

FIG. 30 is a data diagram showing accuracy of variation.

FIG. 31 is a diagram showing an impurity distribution and a variation example thereof according to a conventional simulation system.

FIG. 32 is a diagram showing a process variation and its execution combination by a conventional simulation system.

FIG. 33 is a view showing a part of a search tree of an inference engine by a conventional simulation system.

FIG. 34 is a diagram of optimization of perovskites and Si-based crystals by a new ab initio MD simulator developed by the present inventors.

FIG. 35 is a conceptual configuration diagram of a simulator originally developed by the present inventors.

FIG. 36 is a flowchart of a simulator originally developed by the present inventors.

FIG. 37: P by a simulator developed by the present inventors
The figure of a C4v target example of TO perovskite.

FIG. 38 is a structural diagram of perovskite BaTiO ₃ and PbTiO ₃ .

FIG. 39 is an output diagram of the potential of perovskite PTO by a simulator newly created by the present inventors.

FIG. 40 shows a simulator developed by the present inventors.
Example of making SiO ₂ .

FIG. 41 is a charge output diagram of perovskite PTO by a simulator newly created by the present inventors.

FIG. 42: ab initio / MD newly created by the present inventors
Example of distortion of perovskite PTO by simulator.

FIG. 43: ab initio / MD newly created by the present inventors
Example of oxygen deficiency of perovskite PTO by a simulator.

FIG. 44: ab initio / MD newly created by the present inventors
Example of movement of Pb of perovskite PTO by simulator.

FIG. 45: ab initio / MD newly created by the present inventors
Example of movement of Pb of perovskite PTO by simulator.

FIG. 46 is a part of the simulator data created by the present inventors.

FIG. 47 is a diagram of Sawyer Taw.

FIG. 48 is a diagram of an example of the present invention.

FIG. 49 is a conceptual diagram showing how an in-situ monitor measurement device is added in the clustering tool used in the present invention.

FIG. 50 is a diagram showing how a method of contrasting a limited measurement amount with sufficient information obtained from calculation on the other hand is displayed as a graphic image on a monitor.

FIG. 51 is a diagram showing a temperature sequence as an example of input data to the ab initio molecular dynamics simulation system according to the present invention.

FIG. 52 is a conceptual diagram showing atomic level execution prediction by the ab initio molecular dynamics simulation system according to the present invention.

FIG. 53 is a conceptual diagram visualizing atomic level execution prediction by the abinitio molecular dynamics simulation system according to the present invention.

FIG. 54 is a view showing a time-dependent change of Si—O—Si angle which is an execution example of SiO ₂ by the abinitio molecular dynamics simulation system according to the present invention.

FIG. 55 is a prediction diagram of changes in measured value of temperature, which is an output example for the input sequence of FIG. 51.

FIG. 56 is a distribution diagram of coordination angles in deposited Si which is an example of atomic level execution by the simulation system according to the present invention.

FIG. 57 is a distribution diagram of Si-Si distance, which is an example of atomic level execution by the simulation system according to the present invention.

FIG. 58 is a distribution diagram of the distance between Si and Si, which is an atomic level execution example by the simulation system according to the present invention.

FIG. 59 is a distribution diagram of Si-Si distance, which is an example of atomic level execution by the simulation system according to the present invention.

FIG. 60 is a distribution diagram of Si coordination angles, which is an example of atomic level execution by the simulation system according to the present invention.

FIG. 61 is a distribution diagram of the distance between Si and Si which is an atomic level execution example by the simulation system according to the present invention.

FIG. 62 is an atomic level execution example by the simulation system according to the present invention. The distribution map of Si-Si distance.

FIG. 63 is a distribution diagram of Si coordination angles, which is an atomic level execution example by the simulation system according to the present invention.

FIG. 64 is a distribution diagram of Si coordination angles, which is an atomic level execution example by the simulation system according to the present invention.

FIG. 65 is a distribution diagram of Si coordination angles, which is an example of atomic level execution by the simulation system according to the present invention.

FIG. 66 is a distribution diagram of the distance between Si and Si which is an atomic level execution example by the simulation system according to the present invention.

FIG. 67 is a distribution diagram of the distance between Si and Si which is an atomic level execution example by the simulation system according to the present invention.

FIG. 68 is a distribution diagram of the distance between Si and Si which is an atomic level execution example by the simulation system according to the present invention.

FIG. 69 is a distribution diagram of the distance between Si and Si which is an atomic level execution example by the simulation system according to the present invention.

FIG. 70 is a diagram showing the whole of another system constructed according to the present invention.

FIG. 71 is a characteristic diagram showing conventional Raman measurement results.

FIG. 72 is a characteristic diagram showing conventional Raman measurement results.

FIG. 73 is a characteristic diagram showing conventional Raman measurement results.

FIG. 74 is a diagram showing a model for explaining a conventional Raman measurement result.

FIG. 75 SiO ₄ as a function of O—Si—O bond angle
FIG. 6 is a characteristic diagram showing a wavelength of normal mode vibration of a tetrahedron and a characteristic diagram showing a decrease in strength caused by stress.

FIG. 76 is a model diagram showing defects in SiO ₂ .

77 is a characteristic diagram showing an execution example of the molecular dynamics simulator according to the present invention. FIG.

FIG. 78 is a characteristic diagram showing an execution example of the molecular dynamics simulator according to the present invention.

FIG. 79 is a characteristic diagram showing an execution example of the molecular dynamics simulator according to the present invention.

FIG. 80 is a diagram showing movement of atoms, which is an execution example of the molecular dynamics simulator according to the present invention.

FIG. 81 is a characteristic diagram showing an X-ray spectrum prediction result that is an execution example of the molecular dynamics simulator according to the present invention.

FIG. 82 is a characteristic diagram showing a stress prediction result which is an execution example of the molecular dynamics simulator according to the present invention.

FIG. 83 is a characteristic diagram showing calculation results by a molecular dynamics simulator newly developed by the present inventors.

FIG. 84 is a diagram showing a consistent calculation including a variation of the molecular dynamics simulator using the variational method according to the present invention.

FIG. 85 is a diagram showing a search method by the molecular dynamics simulator according to the present invention.

FIG. 86 is a characteristic diagram showing how temperature is set by the inventors when forming amorphous SiO ₂ .

FIG. 87: Amorphous SiO in a calculator by the present inventors
₂ Characteristic diagram showing the pressure change during creation.

FIG. 88 is a characteristic diagram showing a volume change when an amorphous SiO ₂ is created by a molecular dynamics simulator created by the present inventors.

FIG. 89 is a characteristic diagram showing a volume change when forming amorphous SiO ₂ by a molecular dynamics simulator newly developed by the present inventors.

FIG. 90 is a characteristic diagram showing a potential energy change when forming amorphous SiO ₂ by a molecular dynamics simulator newly developed by the present inventors.

FIG. 91 shows _{2 in} the calculation cell at the time of forming amorphous SiO ₂ by the molecular dynamics simulator newly developed by the present inventors.
The characteristic view which shows the time change of a side.

FIG. 92 is a characteristic diagram showing changes in energy calculated using a normal Verlet algorithm in a molecular dynamics simulator newly developed by the present inventors.

FIG. 93 is a characteristic diagram showing changes in energy calculated by using the improved Verlet algorithm in the molecular dynamics simulator newly developed by the present inventors.

FIG. 94 is a characteristic diagram showing the time change of the dipole moment calculated using the molecular dynamics simulator newly developed by the present inventors.

FIG. 95 is a characteristic diagram showing a power spectrum further obtained from the dipole moment using the molecular dynamics simulator newly developed by the present inventors.

FIG. 96 is a diagram showing an IR spectrum calculated by applying the molecular dynamics simulator newly developed by the present inventors to α-SiO 2.

FIG. 97 is a view showing a predicted IR spectrum when an external pressure is applied by applying the molecular dynamics simulator newly developed by the present inventors to α-SiO 2.

FIG. 98 is a characteristic diagram showing volume change during formation of amorphous SiO by a molecular dynamics simulator newly developed by the present inventors.

[FIG. 99] A normal Verlet algorithm is used for the molecular dynamics simulator newly developed by the present inventors,
The characteristic view which shows the change of the calculated energy.

FIG. 100 is a diagram showing a calculation procedure performed by the present inventors.

101 is a diagram showing a calculation procedure performed by the present inventors. FIG.

102 is a view showing a sliding surface, decomposed shear stress and the like in Si (100). FIG.

FIG. 103 is a schematic diagram showing a reaction in the vicinity of SiO ₂ and a Si / SiO ₂ interface.

FIG. 104 is a diagram schematically showing sequential calculation of viscoelasticity calculation.

FIG. 105 is a characteristic diagram showing temperature-time dependence of the thickness of an oxide film.

FIG. 106 is a characteristic diagram showing temperature-time dependence of the thickness of an oxide film.

FIG. 107 is a characteristic diagram showing temperature-time dependence of the thickness of an oxide film.

FIG. 108 is a diagram showing principal stress distribution and maximum shear stress in trench oxidation when fitted coefficients are used.

FIG. 109 is a characteristic diagram showing changes in temperature and accumulated stress in a typical thermal process.

FIG. 110 is a calculation prediction diagram showing the distribution of σxx at the trench angle.

FIG. 111 is a calculation prediction diagram showing the distribution of σxx at the trench angle.

FIG. 112 is a calculation prediction diagram showing the distribution of σxx at the trench angle.

FIG. 113 is a calculation prediction diagram showing the distribution of σxx at the trench angle.

FIG. 114 is a calculation prediction diagram showing the distribution of σxx at the trench angle.

FIG. 115 is a calculation prediction diagram showing the distribution of σxx at the trench angle.

FIG. 116 is a diagram showing a shape of an optimized oxidation prediction.

117 is a sectional view showing a step of forming a gate thermal oxide film in Example 7. FIG.

118 is a cross-sectional view showing the step of forming the gate thermal oxide film in Example 7. FIG.

FIG. 119 is a characteristic diagram showing a temperature rising step after forming a gate thermal oxide film in Example 7.

FIG. 120 is a characteristic diagram showing changes in stress due to differences in temperature raising steps.

FIG. 121 is a characteristic diagram showing a change in stress due to a difference in heating process.

FIG. 122 is a diagram showing a calculation example of stress distribution.

123 is a diagram showing the relationship between the oxidation temperature and the compressive stress in the film in Example 7. FIG.

FIG. 124 is a diagram showing calculation results of changes in IR spectrum due to stress application.

FIG. 125 is a diagram showing a relationship between a compressive stress in an oxide film and a wave number shift of a Si—O—Si asymmetric expansion / contraction peak.

FIG. 126 is a diagram showing a relationship between Qbd and the appearance frequency of bonds in which the Si—O atom spacing in the first vicinity is 1.8 angstroms or more.

FIG. 127 is a diagram showing changes in Si—O atom spacing distribution with and without stress application.

FIG. 128 is a diagram showing a relationship between dRAM generation and chip cost.

FIG. 129 is a schematic diagram showing a manufacturing process of a semiconductor device according to a conventional example.

FIG. 130 is a diagram showing a clustering tool according to a conventional example.

Continuation of front page (51) Int.Cl. ⁶ Identification number Reference number within the agency FI Technical indication location H01L 29/78 H01L 29/78 301Z 21/336 (72) Inventor Kikuo Yamabe Komukai, Kawasaki City, Kanagawa Prefecture Toshiba-cho 1 Toshiba Research & Development Center (72) Inventor Haruo Okano Komukai Toshiba-cho 1 Kawasaki-shi, Kanagawa 1 Toshiba Research & Development Center

Claims (18)

[Claims] 1. A method of manufacturing a semiconductor device comprising a plurality of steps, wherein a film formed in a manufacturing apparatus in at least one of the plurality of steps, a film formed on a semiconductor substrate, a semiconductor substrate surface, and a semiconductor. A step of measuring at least one of a plurality of characteristics relating to at least one inside the substrate to obtain actual observation data, and an ab initio molecular dynamics process simulator or an empirical potential-provided molecular dynamics simulator, At least one set value of a plurality of manufacturing process factors in at least one, and the fluctuation of this set value as an input value,
A step of performing variational calculation based on this input value to obtain at least one prediction data of the plurality of characteristics; a step of sequentially comparing the prediction data and the actual observation data with each other in real time; Therefore, when a significant difference is found between at least one set value of the plurality of manufacturing process factors and at least one of the plurality of manufacturing process factors inferred from the actual observation data, And a step of sequentially performing correction processing in real time. 2. The ab initio molecular dynamics process simulator or the molecular dynamics simulator to which an empirical potential is applied is provided for the atoms constituting the semiconductor material, insulating material, or conductive material formed or processed in the plurality of steps. When the exercise reaches equilibrium, set the time to 0, and then set the sampling time until time t,
From the abinitio molecular dynamics process simulator or the molecular dynamics simulator given an empirical potential, which has a function of obtaining a time correlation with polarizability, charge, and position vector of each atom and obtaining an optical spectrum at a predetermined time. The method for manufacturing a semiconductor device according to claim 1, wherein the predicted data to be derived includes vibration behavior derived from the momentary motion of each atom, its natural frequency, and distribution information of the natural frequency. . 3. The ab initio molecular dynamics process simulator or the molecular dynamics simulator to which an empirical potential is applied is provided for the atoms constituting the semiconductor material, insulating material, or conductive material formed or processed in the plurality of steps. The material strength constant of the material is sequentially calculated from the motion, and the function of finding the stress induced in the material using this is provided. When the stress exceeds a specified value, at least one of the manufacturing process factors is The method of manufacturing a semiconductor device according to claim 1, wherein one of the semiconductor devices is subjected to a sequential correction process. 4. The method for manufacturing a semiconductor device according to claim 2, wherein the atoms forming the semiconductor material, the insulating material, or the conductive material are silicon atoms and oxygen atoms. 5. A semiconductor device manufacturing apparatus performing a plurality of steps, a film formed in the manufacturing apparatus in at least one of the plurality of steps, a film formed on a semiconductor substrate, a semiconductor substrate surface, and The means for measuring at least one of a plurality of characteristics relating to at least one inside the semiconductor substrate to obtain actual observation data, and the ab initio molecular dynamics process simulator or the empirical potential-provided molecular dynamics simulator, wherein At least one set value of a plurality of manufacturing process factors in at least one of the above, and fluctuations of this set value are used as input values, and variation calculation is performed based on these input values to predict at least one of the plurality of characteristics. Means for obtaining data, means for sequentially comparing and testing the predicted data and actual observation data in real time, and the comparison test Therefore, when a significant difference is found between at least one set value of the plurality of manufacturing process factors and at least one of the plurality of manufacturing process factors inferred from the actual observation data, An apparatus for manufacturing a semiconductor device, comprising: means for sequentially performing correction processing in real time. 6. An ab initio molecular dynamics process simulator or an empirical potential-added molecular dynamics simulator is used to set at least one set value of a plurality of manufacturing process factors in at least one of a plurality of processes of a semiconductor device manufacturing process. , And a fluctuation of this set value as an input value, a variation calculation is performed based on this input value, and a film formed in a manufacturing apparatus in at least one of the plurality of steps, a film formed on a semiconductor substrate. A step of obtaining at least one prediction data of a plurality of characteristics relating to at least one of the semiconductor substrate surface and the inside of the semiconductor substrate, and the prediction data and at least one of the plurality of characteristics in at least one of the plurality of steps. The actual observation data obtained by the measurement is sequentially compared and verified in real time to diagnose the manufacturing process factors and A simulation method comprising a step of making a determination. 7. A significant difference between at least one set value of the plurality of manufacturing process factors and at least one of the plurality of manufacturing process factors inferred from the actual observation data is recognized by the comparison test. 7. The simulation method according to claim 6, further comprising a step of obtaining a corrected value of the manufacturing process factor in real time. 8. At least one of a plurality of characteristics relating to the inside of the manufacturing apparatus, the surface of the semiconductor substrate, and / or the film formed on the semiconductor substrate in at least one of the plurality of steps of the manufacturing process of the semiconductor device is measured. At least one set value of a plurality of manufacturing process factors in at least one of the plurality of processes by means of obtaining actual observation data and an ab initio molecular dynamics process simulator or a molecular dynamics simulator provided with an empirical potential, and The fluctuation of this set value is used as an input value, variational calculation is performed based on this input value, means for obtaining at least one prediction data of the plurality of characteristics, the prediction data and the actual observation data are sequentially and in real time. The simulator for diagnosing and judging the manufacturing process factor of the semiconductor device, which is provided with a means for comparing and certifying the semiconductor device. 9. An ab initio molecular dynamics process simulator or an empirical potential-based molecular dynamics simulator is used to set at least one set value of a plurality of manufacturing process factors in at least one of a plurality of processes of a semiconductor device manufacturing process. , And a fluctuation of this set value as an input value, a variation calculation is performed based on this input value, and a film formed in a manufacturing apparatus in at least one of the plurality of steps, a film formed on a semiconductor substrate. Means for obtaining at least one predictive data of a plurality of characteristics relating to at least one of the semiconductor substrate surface and the inside of the semiconductor substrate, and the predictive data and at least one of the plurality of characteristics in at least one of the plurality of steps. The actual observation data obtained by the measurement is sequentially compared and verified in real time to diagnose the manufacturing process factors and A simulator which comprises a means for making a determination and diagnoses / determines a manufacturing process factor of a semiconductor device. 10. An ab initio molecular dynamics process simulator or an empirical potential-based molecular dynamics simulator is used to set at least one set value of a plurality of manufacturing process factors in at least one of a plurality of processes of a semiconductor device manufacturing process. , And a fluctuation of this set value as an input value, a variation calculation is performed based on this input value, and a film formed in a manufacturing apparatus in at least one of the plurality of steps, a film formed on a semiconductor substrate. Means for obtaining at least one predictive data of a plurality of characteristics relating to at least one of the semiconductor substrate surface and the inside of the semiconductor substrate, and the predictive data and at least one of the plurality of characteristics in at least one of the plurality of steps. The actual observation data obtained by measurement is sequentially compared and verified in real time to diagnose the manufacturing process factors. And a determination method for diagnosing and determining a semiconductor device manufacturing process factor. 11. A field effect type MI on a semiconductor single crystal substrate.
In forming the S element, the ferroelectric film of the MIS element is made of a single crystal, and in a system including the semiconductor substrate and the single crystal ferroelectric film, abini expressing the target characteristic.
An electric field characterized by introducing an evaluation function based on the tio molecular dynamics theory, thereby designing a single crystal ferroelectric film and a substrate according to an optimum solution of a system under desired main operating conditions. Method of manufacturing effective MIS semiconductor device. 12. A field effect type M on a semiconductor single crystal substrate.
In forming an IS element, an ab initio molecule expressing a target characteristic of the system in a system including the semiconductor substrate and the single crystal ferroelectric film, in which the ferroelectric film of the MIS element is a single crystal A field effect type MIS semiconductor, characterized in that the free energy of the entire system is taken as an evaluation function based on the dynamic theory and the single crystal ferroelectric film is arranged so that the free energy is minimized under the main operating conditions. Device manufacturing method. 13. When forming a field effect MIS element on a Si semiconductor single crystal substrate, the gate insulating oxide film of the MIS element is made of single crystal, and the semiconductor substrate and the single crystal ferroelectric film are formed. In the system including, the free energy of the whole system is taken as an evaluation function based on ab initio molecular dynamics theory expressing the target characteristic of the system, and the free energy is minimized under the main operating condition. A method of manufacturing a field-effect SiMIS semiconductor device, comprising disposing a single crystal ferroelectric film. 14. When forming a field effect MIS device on a Si semiconductor single crystal substrate, the gate insulating oxide film of the MIS device is made of single crystal, and the semiconductor substrate and the single crystal ferroelectric film are formed. In the system including, the free energy of the whole system is taken as an evaluation function based on ab initio molecular dynamics theory expressing the object property of the system, and the above free energy is minimized under the main operating condition. A method for manufacturing a field effect type SiMIS semiconductor device, comprising disposing a crystal ferroelectric film and suppressing the oxygen defect concentration of the single crystal ferroelectric film to 0.01% or less. 15. When forming a field effect MIS device on a Si semiconductor single crystal substrate, the gate insulating oxide film of the MIS device is made of single crystal, and the semiconductor substrate and the single crystal ferroelectric film are formed. In the system including, the free energy of the whole system is taken as an evaluation function based on ab initio molecular dynamics theory that expresses the target property of the system, and the free energy of this system is minimized under the main operating conditions. The single crystal ferroelectric film is arranged, and in particular, local unevenness of energy is set within a deviation value 3σ, and the oxygen defect concentration of the single crystal ferroelectric film is suppressed to 0.01% or less. And a method for manufacturing a field effect SiMIS semiconductor device. 16. A method for manufacturing a semiconductor device comprising a plurality of steps, wherein a film formed in a manufacturing apparatus in at least one of the plurality of steps, a film formed on a semiconductor substrate, a semiconductor substrate surface, and a semiconductor. A step of measuring at least one of a plurality of characteristics relating to at least one of the inside of the substrate to obtain actual observation data, and a molecular dynamics simulator given an empirical potential, and a plurality of steps in at least one of the plurality of steps. At least one set value of the manufacturing process factor, and the fluctuation of this setting as an input value, performing variational calculation based on this input value, obtaining at least one prediction data of the plurality of characteristics, and the prediction data And the actual observation data are sequentially and in real time compared and tested, and by the comparison test, at least one of the plurality of manufacturing process factors is If a significant difference is found between the set value and at least one of the plurality of manufacturing process factors inferred from the actual observation data, the manufacturing process factor is sequentially corrected in real time. A method of manufacturing a semiconductor device, comprising: 17. The molecular dynamics simulator to which the empirical potential is applied has a structure in which the motions of atoms forming a semiconductor material, an insulating material, or a conductive material formed or processed in the plurality of steps reach parallel to each other. And set the time to 0
Set to the sampling time until the time t, and has a function of obtaining a time correlation with the charge and position vector of each atom and obtaining an optical spectrum at a predetermined time. The prediction data derived from the science simulator includes vibration behavior derived from the motion of each atom at each time and time, its natural frequency, and distribution information of the natural frequency. Manufacturing method of semiconductor device. 18. The molecular dynamics simulator to which the empirical potential is applied calculates the motion of atoms constituting a semiconductor material, an insulating material or a conductive material formed or processed in the plurality of steps, Input the position vector of atom and electronegativity of each element, or ionization potential and electron affinity of each element, calculate the change in charge due to vibration and displacement of each atom from these data, and then change this charge The method of manufacturing a semiconductor device according to claim 1, wherein the potential is calculated in consideration of the above.