4. Algorithms and techniques for MCS

The MCS is based on a probabilistic process with a random choice of configurations and samples of the situation of the physical system. The two pedagogical examples most cited in the literature are the integration of a single variable function and Ising's model of spin. In the following subsection, we define the integration of a single variable function. We introduce the Ising model at the end of Section 4.2.2.

#### 4.1 Integration of function of a single variable

Calculation of the definite integral for a function f(x) of a single variable x on domain {a, b} has been proposed (Figure 1):

Let:

$$I = \int\_{a}^{b} f(\mathbf{x})d\mathbf{x} \tag{1}$$

Let xi and yi be real random numbers (i = 1, 2,…, N), and let H be a real number greater than the f(x) for x belonging to the domain {a, b} (or x ∈ {a, b}).

Let r1 and r2 be two random numbers belonging to the domain {0, 1} according to a uniform distribution law. Generators (e.g., Ran, RANDOM, RANDUM, or other IMSL mathematical libraries) of random numbers can be used:

Figure 1. The integral of a function f(x).

$$\mathbf{x}\_i = a + r\_1(b - a) \text{ and } y\_i = \mathbf{0} + r\_2(H - \mathbf{0}) \tag{2}$$

b. Choice of an initial configuration that responds to some physical and

How to Use the Monte Carlo Simulation Technique? Application: A Study of the Gas Phase…

thermodynamic properties. The total or internal energy of the system is Ei.

c. Infinitesimal random displacement of a particle (or element of the system) and calculation of the new internal energy of the system Ef. This displacement

is related to the physical magnitudes: time scale and length scale. The physical system tends toward a minimization of the internal energy of

d. If ΔE ≤ 0; the new configuration is retained (favorable) and the different

f. If ε < Pr, accept the move and in any case go back to step (c) for a new choice of an infinitesimal displacement (new configuration). Note that if such a trial move is rejected, the old configuration is again counted in the averaging

Figure 2 shows how to choose between the selected configurations. Let ε be a random number following a uniform law; If ε<sup>1</sup> ≤ Pr the configuration is retained,

Numerical simulation using the MC method is a very important tool for the study of static properties. The basic algorithm is based on probability notions. Understanding of the distribution function and/or interaction potentials is the heart

In equilibrium statistical physics, the system has a certain probability that can be in any states. The probability of being in a state μ with energy H(μ) is given by the

<sup>P</sup><sup>μ</sup> <sup>¼</sup> exp ð Þ �Hð Þ <sup>μ</sup> <sup>=</sup>kBT

<sup>Z</sup> (4)

the system with some fluctuation. Let ΔE=Ef-Ei the fluctuation.

e. If ΔE > 0; a random number ε is chosen such that 0 < ε < 1. Let the probability Pr equal to: Pr = exp. (�ΔE/kBT) (where kB is the Boltzmann

averages can be obtained; go to step (c).

constant and T is the temperature).

DOI: http://dx.doi.org/10.5772/intechopen.88559

with probability Pr.

of the calculation.

Figure 2.

121

Boltzmann distribution P(μ):

and if ε<sup>2</sup> > Pr the configuration is rejected.

4.2.2 Thermodynamic quantities at equilibrium

Configuration choice according to Metropolis scheme.

where xi and yi are random numbers (xi ∈ {a, b} and yi ∈ {0, H}).

The Monte Carlo (MC) method is based on a probabilistic process. Let N be the total number of cases chosen (possible cases). It is necessary to count the number of favorable cases (or the number of points below the curve y = f(x)); let yi ≤ f(xi)). The number of favorable cases is Nfav. When N ➔ ∞, the value I of the integral is [26]:

$$I = \frac{N\_{fav}}{N}(b - a)(H - 0.0) = \frac{N\_{fav}}{N}(b - a)H \tag{3}$$

An example [26] is the calculation of the value π by calculating the integral I on a quarter circle of unit radius (R = 1.0). The pairs of random numbers (xi, yi) satisfying the condition: xi <sup>2</sup> + yi <sup>2</sup> ≤ 1. The function f(x) is equal to ffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>1</sup> � <sup>x</sup><sup>2</sup> <sup>p</sup> .

We take a = 0.0, b = 1.0, and H = 1.0.

For different values of N, we show that the numerical solution tends to π = 4I. Although this integral is simple, it shows the strength and simplicity of the method. The technique can be generalized for the integration of multivariate functions.

We note that integration by the MC method is based on:


## 4.2 Principle of the MCS model

#### 4.2.1 Calculation algorithm (Metropolis algorithm)

For statistical physics problems, the probabilistic choice of configurations is not always deterministic; the favorable and unfavorable cases are not exclusive. According to the Metropolis algorithm [26, 27], the steps of the simulation are:

a. Choice of a simulation cell of adequate shape to the studied phenomena. The size of the simulation cell is related to a scale of length characteristic of the forces and interaction potential of the studied phenomenon. This cell may contain Npc particles (and/or elements).

How to Use the Monte Carlo Simulation Technique? Application: A Study of the Gas Phase… DOI: http://dx.doi.org/10.5772/intechopen.88559


Figure 2 shows how to choose between the selected configurations. Let ε be a random number following a uniform law; If ε<sup>1</sup> ≤ Pr the configuration is retained, and if ε<sup>2</sup> > Pr the configuration is rejected.

Numerical simulation using the MC method is a very important tool for the study of static properties. The basic algorithm is based on probability notions. Understanding of the distribution function and/or interaction potentials is the heart of the calculation.

#### 4.2.2 Thermodynamic quantities at equilibrium

In equilibrium statistical physics, the system has a certain probability that can be in any states. The probability of being in a state μ with energy H(μ) is given by the Boltzmann distribution P(μ):

$$P\_{\mu} = \frac{\exp\left(-H(\mu)/k\_B T\right)}{Z} \tag{4}$$

Figure 2. Configuration choice according to Metropolis scheme.

xi ¼ a þ r1ð Þ b � a and yi ¼ 0 þ r2ð Þ H � 0 (2)

Nfav

<sup>2</sup> ≤ 1. The function f(x) is equal to ffiffiffiffiffiffiffiffiffiffiffiffiffi

<sup>N</sup> ð Þ <sup>b</sup> � <sup>a</sup> <sup>H</sup> (3)

<sup>1</sup> � <sup>x</sup><sup>2</sup> <sup>p</sup> .

where xi and yi are random numbers (xi ∈ {a, b} and yi ∈ {0, H}).

<sup>N</sup> ð Þ <sup>b</sup> � <sup>a</sup> ð Þ¼ <sup>H</sup> � <sup>0</sup>:<sup>0</sup>

<sup>I</sup> <sup>¼</sup> Nfav

We take a = 0.0, b = 1.0, and H = 1.0.

<sup>2</sup> + yi

We note that integration by the MC method is based on:

satisfying the condition: xi

functions.

120

Figure 1.

The integral of a function f(x).

exclusive).

4.2 Principle of the MCS model

4.2.1 Calculation algorithm (Metropolis algorithm)

contain Npc particles (and/or elements).

The Monte Carlo (MC) method is based on a probabilistic process. Let N be the total number of cases chosen (possible cases). It is necessary to count the number of favorable cases (or the number of points below the curve y = f(x)); let yi ≤ f(xi)). The number of favorable cases is Nfav. When N ➔ ∞, the value I of the integral is [26]:

Theory, Application, and Implementation of Monte Carlo Method in Science and Technology

An example [26] is the calculation of the value π by calculating the integral I on a quarter circle of unit radius (R = 1.0). The pairs of random numbers (xi, yi)

For different values of N, we show that the numerical solution tends to π = 4I. Although this integral is simple, it shows the strength and simplicity of the method. The technique can be generalized for the integration of multivariate

• The choice of random configurations according to a uniform distribution law

For statistical physics problems, the probabilistic choice of configurations is not

a. Choice of a simulation cell of adequate shape to the studied phenomena. The size of the simulation cell is related to a scale of length characteristic of the forces and interaction potential of the studied phenomenon. This cell may

always deterministic; the favorable and unfavorable cases are not exclusive. According to the Metropolis algorithm [26, 27], the steps of the simulation are:

• Each configuration chosen is either favorable or unfavorable (the "or" is

where T is the absolute temperature and kB is called Boltzmann's constant. It is conventional to denote the quantity (kBT)�<sup>1</sup> by the symbol β. The normalizing factor Z, or partition function, is given by:

$$Z = \sum\_{\mu} \exp\left(-H(\mu)/k\_B T\right) = \sum\_{\mu} \exp\left(-\beta H(\mu)\right) \tag{5}$$

The Ising model has been studied in one and two dimensions to obtain results of thermal properties, phase transition, and magnetic properties [26–28]. For chosen values of J and/or B, different steps may be taken for the calculations (simulation cell, initialization, configurations, boundary conditions, calculation algorithms). For any configuration, each spin takes the two possible directions. The detail of the

How to Use the Monte Carlo Simulation Technique? Application: A Study of the Gas Phase…

We give a system of N particles (atoms, molecules, ions or particles) placed in a cell of fixed volume, generally of cubic form. The initial positions may, depending on the case, be distributed randomly according to a certain law (uniform or otherwise) or have a given symmetry. In a fluid, a gas, or a plasma, the particles may have random positions in general; in a solid or surface, with a crystal structure, the particles take ordered positions. The choice of random initial positions allows great

At the first step, the particles are given velocities that are generally selected to have a zero total momentum. If the system is in thermodynamic equilibrium, the initial velocities will be randomly chosen according to a Maxwell-Boltzmann law. In the general case, the velocity distribution is according to the problem dealt with. All other phase properties can be initialized to the particles; the main thing is the

The particles interact with each other according to chosen interaction potentials. Since the interaction potentials are specific for each "numerical experiment," the main part of the work consists in calculating the interaction energies for each

The choice of interaction potentials is directly related to the mathematical formulation of the problem according to the state of the medium: fluid, gas, plasma, or solid. It can be Lennard-Jones potential, Coulomb potential, Debye potential, Morse potential, Stillinger-Weber potential, Born-Mayer potential, Moliere potential, or

In general, two main boundary conditions are used: periodic boundary condi-

To minimize the surface effect, periodic boundary conditions (PBC) [30] are invariably imposed. The simulation cell is reproduced throughout the space to form an infinite mesh. We can simulate the properties of an infinite system. The particles that we follow are in the central cell; if a particle crosses a wall with a certain velocity, its image returns with the same velocity by the opposite wall. Under these conditions, the number of particles in the central cell, and consequently the density, is constant. These conditions also allow the conservation of the energy and the momentum of the system and do not introduce periodic effects (because of the

According to the hypotheses and according to the geometry of the problem, other boundary conditions are proposed [26]. For example, in order to model thin films, the simulation cells are longitudinal and parallel to the film; one uses PBC in

calculation procedure is not the purpose of this chapter.

freedom on the choice of the number of particles in the cell.

conservation of the total quantities of the system.

tions (PBC) and minimum image convention (MIC) [29].

4.2.3.2 Potentials of interaction

proposed configuration.

4.2.3.3 Boundary conditions

interaction between particles).

others.

123

4.2.3 MCS module designs

4.2.3.1 Simulation cell and initialization

DOI: http://dx.doi.org/10.5772/intechopen.88559

The average of a quantity Q fora system in equilibrium is:

$$ = \sum\_{\mu} Q\_{\mu} P\_{\mu} = \frac{1}{Z} \sum\_{\mu} Q\_{\mu} \exp\left(-\beta H(\mu)\right) \tag{6}$$

The internal energy U, is given by:

$$U = \frac{1}{Z} \sum\_{\mu} H(\mu) \exp\left(-\beta H(\mu)\right) \tag{7}$$

which can be written in terms of a derivative of the partition function:

$$U = \frac{1}{Z} \frac{\partial Z}{\partial \beta} = -\frac{\partial \log Z}{\partial \beta} \tag{8}$$

From thermodynamics we have expressions for the specific heat C, the entropy S, and the Helmholtz free energy F:

$$\mathbf{C} = \frac{\partial U}{\partial T} = -k\_B \beta^2 \frac{\partial U}{\partial \beta} = -k\_B \beta^2 \frac{\partial^2 \log Z}{\partial \beta^2} \tag{9}$$

or

$$\mathbf{C} = T \frac{\partial \mathbf{S}}{\partial T} = -\beta \frac{\partial \mathbf{S}}{\partial \beta} \tag{10}$$

and

$$S = -k\_B \beta \frac{\partial \log Z}{\partial \beta} + k\_B \log Z \tag{11}$$

and

$$F = U - T\mathbf{S} = -k\_B \log Z \tag{12}$$

We can calculate other parameters affecting the system.

The Monte Carlo method is an excellent technique for estimating probabilities, and we can take advantage of this property in evaluating the results. The simplest and most popular model of a system of interacting variables in statistical physics is the Ising model. It consists of spins σ<sup>i</sup> which are confined to the sites of a lattice and which may have only the values (+1) and (�1). These spins interact with their nearest neighbors on the lattice with interaction constant J; they can interact with an external magnetic field B coupling to the spins. The Hamiltonian H for this model is [26]:

$$H = -J\sum\_{i,j} \sigma\_i \sigma\_j - B\sum\_i \sigma\_i \tag{13}$$

How to Use the Monte Carlo Simulation Technique? Application: A Study of the Gas Phase… DOI: http://dx.doi.org/10.5772/intechopen.88559

The Ising model has been studied in one and two dimensions to obtain results of thermal properties, phase transition, and magnetic properties [26–28]. For chosen values of J and/or B, different steps may be taken for the calculations (simulation cell, initialization, configurations, boundary conditions, calculation algorithms). For any configuration, each spin takes the two possible directions. The detail of the calculation procedure is not the purpose of this chapter.

#### 4.2.3 MCS module designs

where T is the absolute temperature and kB is called Boltzmann's constant. It is conventional to denote the quantity (kBT)�<sup>1</sup> by the symbol β. The normalizing

μ

exp ð Þ �βHð Þ μ (5)

Q<sup>μ</sup> exp ð Þ �βHð Þ μ (6)

<sup>∂</sup><sup>β</sup> (8)

<sup>∂</sup>β<sup>2</sup> (9)

<sup>∂</sup><sup>β</sup> (10)

σ<sup>i</sup> (13)

<sup>∂</sup><sup>β</sup> <sup>þ</sup> kB log <sup>Z</sup> (11)

F ¼ U � TS ¼ �kB log Z (12)

Hð Þ μ exp ð Þ �βHð Þ μ (7)

exp ð Þ¼ �Hð Þ <sup>μ</sup> <sup>=</sup>kBT <sup>X</sup>

Theory, Application, and Implementation of Monte Carlo Method in Science and Technology

<sup>Q</sup>μP<sup>μ</sup> <sup>¼</sup> <sup>1</sup> Z X μ

which can be written in terms of a derivative of the partition function:

<sup>∂</sup><sup>β</sup> ¼ � <sup>∂</sup> log <sup>Z</sup>

<sup>∂</sup><sup>T</sup> ¼ �<sup>β</sup> <sup>∂</sup><sup>S</sup>

<sup>∂</sup><sup>β</sup> ¼ �kBβ<sup>2</sup> <sup>∂</sup><sup>2</sup> log <sup>Z</sup>

From thermodynamics we have expressions for the specific heat C, the entropy

<sup>U</sup> <sup>¼</sup> <sup>1</sup> Z ∂Z

<sup>∂</sup><sup>T</sup> ¼ �kBβ<sup>2</sup> <sup>∂</sup><sup>U</sup>

<sup>C</sup> <sup>¼</sup> <sup>T</sup> <sup>∂</sup><sup>S</sup>

<sup>S</sup> ¼ �kB<sup>β</sup> <sup>∂</sup> log <sup>Z</sup>

The Monte Carlo method is an excellent technique for estimating probabilities, and we can take advantage of this property in evaluating the results. The simplest and most popular model of a system of interacting variables in statistical physics is the Ising model. It consists of spins σ<sup>i</sup> which are confined to the sites of a lattice and which may have only the values (+1) and (�1). These spins interact with their nearest neighbors on the lattice with interaction constant J; they can interact with an external magnetic field B coupling to the spins. The Hamiltonian H for this model is [26]:

We can calculate other parameters affecting the system.

H ¼ �J

X i,j

σiσ<sup>j</sup> � B

X i

The average of a quantity Q fora system in equilibrium is:

μ

<sup>U</sup> <sup>¼</sup> <sup>1</sup> Z X μ

<sup>C</sup> <sup>¼</sup> <sup>∂</sup><sup>U</sup>

factor Z, or partition function, is given by:

<sup>Z</sup> <sup>¼</sup> <sup>X</sup> μ

The internal energy U, is given by:

S, and the Helmholtz free energy F:

or

and

and

122

<sup>&</sup>lt; <sup>Q</sup><sup>&</sup>gt; <sup>¼</sup> <sup>X</sup>

#### 4.2.3.1 Simulation cell and initialization

We give a system of N particles (atoms, molecules, ions or particles) placed in a cell of fixed volume, generally of cubic form. The initial positions may, depending on the case, be distributed randomly according to a certain law (uniform or otherwise) or have a given symmetry. In a fluid, a gas, or a plasma, the particles may have random positions in general; in a solid or surface, with a crystal structure, the particles take ordered positions. The choice of random initial positions allows great freedom on the choice of the number of particles in the cell.

At the first step, the particles are given velocities that are generally selected to have a zero total momentum. If the system is in thermodynamic equilibrium, the initial velocities will be randomly chosen according to a Maxwell-Boltzmann law. In the general case, the velocity distribution is according to the problem dealt with. All other phase properties can be initialized to the particles; the main thing is the conservation of the total quantities of the system.

#### 4.2.3.2 Potentials of interaction

The particles interact with each other according to chosen interaction potentials. Since the interaction potentials are specific for each "numerical experiment," the main part of the work consists in calculating the interaction energies for each proposed configuration.

The choice of interaction potentials is directly related to the mathematical formulation of the problem according to the state of the medium: fluid, gas, plasma, or solid. It can be Lennard-Jones potential, Coulomb potential, Debye potential, Morse potential, Stillinger-Weber potential, Born-Mayer potential, Moliere potential, or others.

#### 4.2.3.3 Boundary conditions

In general, two main boundary conditions are used: periodic boundary conditions (PBC) and minimum image convention (MIC) [29].

To minimize the surface effect, periodic boundary conditions (PBC) [30] are invariably imposed. The simulation cell is reproduced throughout the space to form an infinite mesh. We can simulate the properties of an infinite system. The particles that we follow are in the central cell; if a particle crosses a wall with a certain velocity, its image returns with the same velocity by the opposite wall. Under these conditions, the number of particles in the central cell, and consequently the density, is constant. These conditions also allow the conservation of the energy and the momentum of the system and do not introduce periodic effects (because of the interaction between particles).

According to the hypotheses and according to the geometry of the problem, other boundary conditions are proposed [26]. For example, in order to model thin films, the simulation cells are longitudinal and parallel to the film; one uses PBC in the directions parallel to the film. In the direction normal to the film, free edge boundary conditions can be used. In such cases, it may be appropriate to also include surface fields and surface interactions. In this way, one can study phenomena such as wetting, interface localization-delocalization transitions, surfaceinduced ordering and disordering, etc.

The core of the program includes calculating the potential energies of particle configuration and particle collisions. The interactions and collisions between particles can be elastic or inelastic; they can be binary or collective. For computation, the interaction energy of a particle with its neighbors is carried out by refocusing a base cell on the particle. This particle only interacts with particles in this region. This is called the "minimal image convention" (MIC) [1].

### 4.2.3.4 Sampling of random data

Generally, a RANDOM generator of real random numbers ri belonging to the domain {0, 1} (or ri ∈ {0, 1} is available. This distribution law is uniform.

To have a real random number xi belonging to the domain {a, b} (or xi ∈ f g a; b ) according to a law of uniform distribution, we have:

$$x\_i = a + r\_i(b - a) \tag{14}$$

corresponds to a random value xran of the domain {xj-1, xj

an irregular mesh Δxi, or with tabular data f(xi) with i = 1,…, m.

This technique can be generalized for a nonuniform distribution law f(x) with

How to Use the Monte Carlo Simulation Technique? Application: A Study of the Gas Phase…

The technique can be generalized, too, for a discrete distribution law f(i) with

In the literature, the reader can find simple algorithms for the choice of random

It is necessary to find some parameters allowing the control of the smooth course of the evolution of the system. We must look for the constants of movement. For example for an isolated system, we have the conservation of the total energy and the

By using the numerical simulation, it is possible to calculate many spatiotemporal quantities F(r,t). These quantities can be positions, speeds, kinetic moments, particle energies, concentrations, transport coefficients, etc. It would then be pos-

For the calculation of the averages, one can note the quantities on the space, on the time or on both. The histogram methods can be used. Static or dynamic distribution functions and spatial or temporal correlation functions can be calculated. It should be noted that the SMC is much more adequate for static properties because

Any calculated function or parameter F(r,t) can be used for another application

In the MCS model discussed extensively in this chapter, it's more about collisions

between particles. It's particle-particle MCS or PP-MCS. In many problems of physics, the general idea is the same, but the applications and proposed models are

Other MCS models, named particle-in-cell MCS (PIC-MCS), are based on particle-cell interactions. In these last models, we also use a probabilistic choice of configurations and small variations in the state of the system (following the

Metropolis algorithm); the interaction is between the particle with a cell, a mesh, or a drop. The parameters and variables of the cell, although local and instantaneous, are macroscopic. These parameters and variables can be thermodynamic, fluid, or electromagnetic. An example of the model based on PIC-MCS is described by Mattei et al. [31] for simulation of electromagnetic particle-in-cell collision in inductively coupled plasmas. Several works can be found in the literature on this same line of

For statistical physics problem solving (such as thin film deposition problems), MCS models use experimental, numerical, or theoretical data from other methods and models. Models can be improved to hybrid models. In the hybrid models, connections between two modules can be realized. The first module is MCS; the second module is fluid, electromagnetic, or other. An example of a three-module hybrid model is presented by Mao and Bogaerts [33] to study gas mixtures in PECVD system. The three modules are MCS, fluid, and electromagnetic. The first

work. Other MCS models using particles may be considered. [32].

formula (or the law) of nonuniform distribution f(x).

DOI: http://dx.doi.org/10.5772/intechopen.88559

numbers of some simple functions (Gaussian, etc.).

4.2.3.5 Control of the evolution of the physical system

sible to calculate all other quantities related to F(r,t).

4.3 Other large methods of Monte Carlo simulation

of the probabilistic choice of configurations.

in another calculation program.

numerous.

125

i = 1,…, m.

quantity of matter.

4.2.3.6 Statistical calculations

}; this number satisfies the

To have a real random number xi belonging to the domain {a, b} (or xi ∈ {a,b}) according to a formula (or law) of nonuniform distribution f(x), a histogram technique is used. Let Nm be the number of intervals. If the mesh is regular (Figure 3):

$$
\Delta \mathbf{x} = (\mathbf{b} - \mathbf{a}) / \mathbf{N}\_{\mathbf{m}} \tag{15}
$$

We define:

$$\mathbf{f}\_{\mathbf{i}} = \mathbf{f}(\mathbf{x}\_{\mathbf{i}}) \text{ for } \mathbf{i} = \mathbf{0}, \dots, \mathbf{m} \text{ and } \mathbf{:} \ \mathbf{x}\_{\mathbf{i}} = \mathbf{a} + \mathbf{i}. \Delta \mathbf{x} \tag{16}$$

We define the sequence:

$$\mathbf{S}\_{\mathbf{0}} = \mathbf{0} \text{ and } \mathbf{:S}\_{\mathbf{i}+\mathbf{1}} = \mathbf{S}\_{\mathbf{i}} + \Delta \mathbf{x} (\mathbf{f}(\mathbf{x}\_{\mathbf{i}}) + \mathbf{f}(\mathbf{x}\_{\mathbf{i}-1}))/2 \tag{17}$$

and the sequence:

$$\mathbf{r}\mathbf{x}\_0 = \mathbf{0} \text{ et } \mathbf{r}\mathbf{x}\_i = \mathbf{S}\_i/\mathbf{S}\_{\mathbf{m}} \tag{18}$$

Hence each real random number ri belongs to the domain {0, 1} (where ri ∈ {0, 1}) (according to the uniform law); this number belongs to the domain {rxj-1, rxj }. It

Figure 3. Random number selection according to f (x) distribution.

How to Use the Monte Carlo Simulation Technique? Application: A Study of the Gas Phase… DOI: http://dx.doi.org/10.5772/intechopen.88559

corresponds to a random value xran of the domain {xj-1, xj }; this number satisfies the formula (or the law) of nonuniform distribution f(x).

This technique can be generalized for a nonuniform distribution law f(x) with an irregular mesh Δxi, or with tabular data f(xi) with i = 1,…, m.

The technique can be generalized, too, for a discrete distribution law f(i) with i = 1,…, m.

In the literature, the reader can find simple algorithms for the choice of random numbers of some simple functions (Gaussian, etc.).

#### 4.2.3.5 Control of the evolution of the physical system

It is necessary to find some parameters allowing the control of the smooth course of the evolution of the system. We must look for the constants of movement. For example for an isolated system, we have the conservation of the total energy and the quantity of matter.

#### 4.2.3.6 Statistical calculations

the directions parallel to the film. In the direction normal to the film, free edge boundary conditions can be used. In such cases, it may be appropriate to also include surface fields and surface interactions. In this way, one can study phenomena such as wetting, interface localization-delocalization transitions, surface-

Theory, Application, and Implementation of Monte Carlo Method in Science and Technology

The core of the program includes calculating the potential energies of particle configuration and particle collisions. The interactions and collisions between particles can be elastic or inelastic; they can be binary or collective. For computation, the interaction energy of a particle with its neighbors is carried out by refocusing a base cell on the particle. This particle only interacts with particles in this region. This is

Generally, a RANDOM generator of real random numbers ri belonging to the

To have a real random number xi belonging to the domain {a, b} (or xi ∈ f g a; b )

To have a real random number xi belonging to the domain {a, b} (or xi ∈ {a,b})

according to a formula (or law) of nonuniform distribution f(x), a histogram technique is used. Let Nm be the number of intervals. If the mesh is regular

xi ¼ a þ rið Þ b � a (14)

Δx ¼ ð Þ b � a =Nm (15)

rx0 ¼ 0 et rxi ¼ Si=Sm (18)

}. It

fi ¼ f xð Þ<sup>i</sup> for i ¼ 0, …:, m and : xi ¼ a þ i: Δx (16)

S0 ¼ 0 and : Siþ<sup>1</sup> ¼ Si þ Δxð Þ f xð Þþ<sup>i</sup> f xð Þ <sup>i</sup>�<sup>1</sup> =2 (17)

Hence each real random number ri belongs to the domain {0, 1} (where ri ∈ {0, 1})

(according to the uniform law); this number belongs to the domain {rxj-1, rxj

domain {0, 1} (or ri ∈ {0, 1} is available. This distribution law is uniform.

induced ordering and disordering, etc.

4.2.3.4 Sampling of random data

(Figure 3):

Figure 3.

124

We define:

We define the sequence:

Random number selection according to f (x) distribution.

and the sequence:

called the "minimal image convention" (MIC) [1].

according to a law of uniform distribution, we have:

By using the numerical simulation, it is possible to calculate many spatiotemporal quantities F(r,t). These quantities can be positions, speeds, kinetic moments, particle energies, concentrations, transport coefficients, etc. It would then be possible to calculate all other quantities related to F(r,t).

For the calculation of the averages, one can note the quantities on the space, on the time or on both. The histogram methods can be used. Static or dynamic distribution functions and spatial or temporal correlation functions can be calculated. It should be noted that the SMC is much more adequate for static properties because of the probabilistic choice of configurations.

Any calculated function or parameter F(r,t) can be used for another application in another calculation program.

#### 4.3 Other large methods of Monte Carlo simulation

In the MCS model discussed extensively in this chapter, it's more about collisions between particles. It's particle-particle MCS or PP-MCS. In many problems of physics, the general idea is the same, but the applications and proposed models are numerous.

Other MCS models, named particle-in-cell MCS (PIC-MCS), are based on particle-cell interactions. In these last models, we also use a probabilistic choice of configurations and small variations in the state of the system (following the Metropolis algorithm); the interaction is between the particle with a cell, a mesh, or a drop. The parameters and variables of the cell, although local and instantaneous, are macroscopic. These parameters and variables can be thermodynamic, fluid, or electromagnetic. An example of the model based on PIC-MCS is described by Mattei et al. [31] for simulation of electromagnetic particle-in-cell collision in inductively coupled plasmas. Several works can be found in the literature on this same line of work. Other MCS models using particles may be considered. [32].

For statistical physics problem solving (such as thin film deposition problems), MCS models use experimental, numerical, or theoretical data from other methods and models. Models can be improved to hybrid models. In the hybrid models, connections between two modules can be realized. The first module is MCS; the second module is fluid, electromagnetic, or other. An example of a three-module hybrid model is presented by Mao and Bogaerts [33] to study gas mixtures in PECVD system. The three modules are MCS, fluid, and electromagnetic. The first

value of densities of species. It provides information on SiH4 dissociation and on the

How to Use the Monte Carlo Simulation Technique? Application: A Study of the Gas Phase…

take into account electrons and eight neutral species (SiH4, SiH3, SiH2, H, H2, Si2H6, Si2H5, SiH). Reactions taken into account include seven electron-neutral and 14 neutral-neutral reactions. Table 1 shows the 21 reactions and rate constants Kreac. At low temperature, the neutrals interact occasionally with each other and move under the effect of thermal agitation; their velocity distribution function is Maxwell-Boltzmann distribution. Electrons have the mean velocity with kinetic energy Te.

reac

!<sup>K</sup>prod reac c <sup>0</sup> � A þ d<sup>0</sup>

R7 Si2H6 + e!SiH4 + SiH2 +e K7 = 1.1 � <sup>10</sup>10� (1.(1./(1. + (0.63 � P)))) [43] R8 SiH4 + H!SiH3 + H2 K8 = 2.8 � <sup>10</sup>�<sup>11</sup> � exp.(�1250/T) [44] R9 SiH4 + SiH2!Si2H6 K9 = 1.1 � <sup>10</sup><sup>10</sup> � (1.�(1./(1. + (0.63 � P)))) [43]

R10 SiH3 + SiH3!SiH4 + SiH2 K10 = 0.45 � 1.5 � <sup>10</sup>�<sup>10</sup> [44] R11 SiH4 + Si2H5!SiH3 + Si2H6 K11 = 5 � <sup>10</sup>�<sup>13</sup> [42] R12 SiH3 + H!SiH2 + H2 K12 = 2 � <sup>10</sup>�<sup>11</sup> [44]

R13 SiH3 + Si2H6!SiH4 + Si2H5 K13 = 4 � <sup>10</sup>�<sup>10</sup> � exp. (�2500/T) [44]

R15 Si2H6 + H!Si2H5 + H2 K15 = 0.66 � 2.4 � <sup>10</sup>�<sup>10</sup> � exp. (�1250/T) [43] R16 Si2H6 + H!SiH4 + SiH3 K16 = 0.34 � 2.4 � <sup>10</sup>�<sup>10</sup> � exp. (�1250/T) [44]

R19 SiH2 + H2!SiH4 K19 = 3 � <sup>10</sup>�<sup>12</sup> � (1. + (1./1. + (0.03 � P))) [43]

R21 SiH4 + SiH!Si2H5 K21 = (1.�(1./(1. + (0.33 � P)))) � (6.9 � <sup>10</sup>�10) [43]

R14 SiH2 + H!SiH + H2 k14 = 2 � <sup>10</sup>�<sup>11</sup> [44]

R17 SiH + H2!SiH3 K17 = 2 � <sup>10</sup>�<sup>12</sup> [43] R18 SiH2 + SiH3!Si2H5 K18 = 3.77 � <sup>10</sup>�<sup>13</sup> [43]

R20 2SiH3!Si2H6 K20 = 0.1 � 1.5 � <sup>10</sup>�<sup>10</sup> [43]

List of gas phase reactions and corresponding rate constants [24].

The plasma in the PECVD reactor is weakly ionized. For our study, the mixture gas contains 22% of SiH4 and 78% of H2; the pressure is 100 mtorr, the temperature of the gas ranges from 373 to 723 K, the electron temperature is about 2.5 eV, and

. The process is considered to be stationary. We

reac be the rate constants of the consumption and the production

c � C þ d � D

� D<sup>0</sup>

/s)

production of SiH3, H, SiH2, and Si2H6 and other important parameters.

of species A. The chemical reaction for the consumption of A is as:

And chemical reaction for the production of A is as:

a<sup>0</sup> � A<sup>0</sup> þ b<sup>0</sup> � B<sup>0</sup>

Symbol Reactions Kreac (cm3

R1 SiH4 + e!SiH3 + H+e k1 = 3 � <sup>10</sup>�<sup>11</sup> [42] R2 SiH4 + e!SiH2 + 2H + e K2 = 1.5 � <sup>10</sup>�<sup>10</sup> [42] R3 SiH4 + e!SiH + H + H2 + e K3 = 9.34 � <sup>10</sup>�<sup>12</sup> [42] R4 SiH4 + e!SiH2 + H2 + e K4 = 7.19 � <sup>10</sup>�<sup>12</sup> [42] R5 H2 + e!2H + e K5 = 4.49 � <sup>10</sup>�<sup>12</sup> [42] R6 Si2H6 + e!SiH3 + SiH2 +H+e K6 = 3.72 � <sup>10</sup>�<sup>10</sup> [42]

<sup>a</sup> � <sup>A</sup> <sup>þ</sup> <sup>b</sup> � <sup>B</sup>!<sup>K</sup>cons

the electron density is 3. 108 cm�<sup>3</sup>

DOI: http://dx.doi.org/10.5772/intechopen.88559

reac and Kprod

Let Kcons

Table 1.

127

### Figure 4.

Schematic of a hybrid model of three modules used to study gas mixtures in the PECVD [33].

module EM calculates the electromagnetic fields by solving Maxwell equations. These fields are used as inputs in the module MCS, where the electron density, electron temperature, electron energy distribution function, and electron impact reaction rates can be computed with a Monte Carlo procedure. Subsequently, the module fluid calculates densities and fluxes of the various plasma species (i.e., heavy particles and electrons) with continuity equations and the electrostatic field with Poisson's equation. This electrostatic field is used as input again in the EM. This cycle is iterated until convergence. The schematic of the hybrid model is given in Figure 4.

To solve statistical physics problems with evolutions as a function of time, kinetic models of MCS (kMCS) are used. Using kMCS, Battaile and Srolovitz [17] described kinetic phenomena of the diffusive motion of a single interstitial atom in a close-packed metal crystal. The motion of the interstitial atom is usually limited to two types: vibration of the atom around the center of the interstitial hole in which it resides and hops to nearest-neighbor interstitial sites. The atom can hop into any of the nearest-neighbor interstitial sites; it executes a random walk. In an MC simulation of this diffusion process, the new position of the interstitial atom is chosen at random from a list of the adjacent interstitial sites.

Other CVD and PECVD works on MCS are presented in Ref.s [15, 34–38]. They show how MCS methods can study properties of gas mixtures and properties of the growth of thin films.
