Simulating Chemical Vapor Deposition

This tutorial

  • explains the modeling concepts behind simulating chemical vapor deposition
  • provides a step-by-step guide on how to set up such a simulations from the GUI

Modeling Chemical Vapor Deposition

From a modeling perspective, chemical vapor deposition (CVD) or atomic layer deposition (ALD) is similar to other processes involving a surface interacting with an atmosphere above it. Further examples of surface-atmosphere interactions include physical vapor deposition (PVD), solid-gas phase heterogeneous catalysis, as well as plasma-surface reactions. As such all of these processes can in principle be studied with setups similar to the CVD example application introduced below.


Simulating surface-atmosphere interactions typically involves quite large computations due to the following reasons.

Model Size
The employed periodic surface models have to feature large unit cells to account for the amorphous structural motifs that result from the deposition of nanometer sized films. In the case of more energetic particles (plasma, sputtering, etc.) impacting the surface, the slab model of the latter needs to be have a certain thickness to account for particle penetration phenomena.
Time Scale
To model reactions on a surface the time resolution of the simulation needs to be short enough to account for the highest vibrational modes of the system. This time scale is often orders of magnitude smaller than the average interval between successive particle impacts around any given site of the model. When studying atmospheres with low pressures such as those occuring in typical high vacuum surface science experiments, this discrepancy in time scales increases further. Furthermore, surface reactions may involve several intermediate or diffusion steps and require therefore significant simulation times too.
Even with the accelerated dynamics simulations over millions of steps presented in this tutorial, the effective pressures achieved during the simulation are much higher than in reality. While being a compromise, this approach can be justified for cases where the chemical evolution at any given site occurs fast enough to not be affected by successive particle impacts.
The local temperature of the slab model might be raised significantly by either a too strong flux of particles, or even by some acceleration methods. The simulation needs thus to be monitored for temperature spikes and adjusted in order to avoid unrealistic side reactions.
The quality with which a model is described during its simulation is crucial to study a realistic evolution of the system. While the model description of such long simulations is essentially limited to empirical force fields or machine learning potentials, it is advised to carefully parametrize such methods for relevant reaction steps prior to begin any production simulations. While the example discussed below simulates realistic chemical processes with an existing ReaxFF parametrization, this cannot be assumed in the general case. Especially, relative discrepancies in reaction barrier heights can affect the chemical processes occurring during the simulation and significantly alter the evolution of the system. It is therefore the responsibility of the user to use adequately validated model descriptions for the simulations in interest. The params tutorials and documentation provide further information on how to obtain such parameterizations.


The model used in the CVD simulations is sketched below and features the following components


Fig. 26 Sketch of the CVD simulation model: particles are inserted in a box (orange) and translated close to their first impact point, where they are actually placed into the simulation box. From there, a particle can either be adsorbed on the surface and undergo further reactions there (blue), or either be reflected (red) or penetrate through the surface (green). In either of the latter cases, it will leave the safe box (dashed line) and be removed.

Mechanical Properties
Our model mimics a surface with an underlying bulk, and should therefore have a similar mechanical resistance to particle impacts. This is achieved by fixing all atoms in a mechanical support layer at the bottom of the slab model.
Thermal Properties
The model should also approach the thermal properties of the underlying bulk material while still allowing for local hotspots around particle impact sites. This is implemented by adding a thermostat module affecting only atoms in another layer deep below the surface, directly above the mechanical support.
Molecule Gun
Particles are (formally) inserted in a random position within a zone located several nanometers above the surface and with a velocity vector pointing downwards.
Particle Translation
To save simulation steps during the particle flight time, the particle is translated close to the first impact point in its flight path. It is then actually placed into the simulation box near this impact point.
Accelerated Dynamics
In regular intervals timed to start directly after each intial particle contact with the surface, a sequence of accelerated dynamics steps using the force-bias Monte Carlo method commences.
Molecule Sink
If a particle moves away from the surface, it will likely not interact with any other parts of the model either. A safe-box is defined encompassing a large region around the entire surface model. If any stray particle leaves the safe box it is automatically removed by the molecule sink module.

Simulation Setup and Settings

Oxidized Ge(001) Surface

The slab model was obtained from a crystal Ge structure, whose lattice parameters were optimized with the ReaxFF parametrization CHOFeAlNiCuSCrSiGe.ff used in all subsequent simulations. After converting the crystal to a surface slab model, the oxygen centers were added to the surface and relaxed in the following. We refer to the crystals & surfaces tutorial for details on building such slab models and provide the initial structure of the partially oxidized Ge(001) surface as an xyz-file for download.

2. In AMSinput select File Import Coordinates and select the downloaded file


With the initial model imported, we begin by setting up the mechanical support and thermostat layers in the model:

1. Switch to a lateral view of the model (press Ctrl + 1 or Ctrl + 2 )
2. Select the lowest (001)-layer of atoms (see Figure below)
3. Change to the Model Regions panel, click AddButton and name the new region as fixed
4. Repeat the above steps for the second lowest (001)-layer and name the region thermo
5. Change to Model Geometry Constraints and PES Scan, select one atom from the fixed region
6. Click on AddButton fixed (fixed position)

Fig. 27 Mechanical support and thermostat layer as regions used in the simulations


We can now set up a short simulation to equilibrate the system to the temperature used in the following calculations.

1. Switch to engine ReaxFFPanel and select Force Field: CHOFeAlNiCuSCrSiGe.ff
2. Set Task: Molecular Dynamics and on MoreBtn
3. Set Number of steps: 100000, Sample Frequency: 100, Initial Temperature 500 K (optional: Checkpoint Frequency 100000000)

Fig. 28 Molecular dynamics settings

4. Change to the thermostat details by clicking on the MoreBtn next to it
5. Click on AddButton to add a thermostat instance
6. Set Thermostat: NHC, Temperature(s) 500 K, Damping constant: 250 fs, Atoms in region: thermo

Fig. 29 Thermostat settings

While the above settings suffice to equilibrate the system in principle, we would also like to include the effect of the accelerated dynamics phases used in the main simulation. Indeed, the employed force bias Monte Carlo method is known to accelerate the evolution of a system by violating the reversibility of larger Monte Carlo steps. For this reason, we activate the fbMC module already for the equilibration run:

7. Change to Model Force Bias MC (fbMC)
8. Click on AddButton to add a new fbMC instance
9. Set Frequency: 12000, Number of steps: 3000, Step length: 0.06 Å, Temperature: 450 K, Enable molecular moves Yes, Translation step length: 0.05 Å, Rotation step angle: 0.05 Radian

Fig. 30 Force bias Monte Carlo settings


We deliberately use a slightly lower temperature during the fbMC phase than for the thermostat, 450 vs 500 K, to avoid too high temperature spikes during the simulation. This choice is justified as we are interested in the general chemical evolution of the system. In addition, because the impinging particle flux adds energy to the system, keeping the temperature stable is deemed more important in this situation than a precise thermodynamic sampling. In other applications, the effect of the fbMC temperature and the fbMC step length on the temperature actually observed need to be validated by the user.

10. Save the file under File Save and run the calculation with File Run


Using the equilibrated positions and velocities allows now to set up the production simulation. We load

after the above equilibration run is complete

1. (Re)open the equilibration calculation in AMSinput
2. In the popup window, enable Use MD velocities from results for AMS MD restart
3. Click on Yes, new job
4. File Save as… and save the input under a new name


Make sure that the original job and the ams.rkf result file is not overwritten

For convenience, we provide the input file for the equilibration run, Ge200O25_equilibration.ams, as well as the resulting Ge200O25.rkf file containing the equilibrated Ge(001) surface structure and atomic velocities for download.

1. Download the files Ge200O25.rkf and Ge200O25_equilibration.ams
2. Open AMSinput
3. File Import Coordinates… and select Ge200O25.rkf


In newer versions of AMS, the regions are imported from the .rkf file. If the regions fixed and thermo are not present after importing the coordinates, either update AMS or repeat the steps in the above Regions subsection.

4. File Import Settings… and select Ge200O25_equilibration.ams
5. Switch to Models MD
6. Set Initial velocities: From File
7. File: Ge200O25.rkf

After either of these steps, we can start setting up the final simulation run. We begin with adapting the settings for MD and fbMC

1. Switch to Models MD
2. Set Number of steps: 15000000
3. Switch to Models fbMC
4. Set Start Step: 7000

Note that the above Start Step settings are used to align the start of the fbMC intervals with the molecule gun shots.

We proceed by loading the structures for the O2 and AlH3 projectiles. These are provided in the files and for download.

1. Download the files and
2. Switch to Models Regions
3. File Import Coordinates… and select
4. Add the imported atoms to a new region and name it O2
5. Repeat the above steps for the file and name the corresponding region AlH3

We can now use these imported regions as O2 and AlH3 projectiles in the molecule gun module.

1. Switch to Models Molecule Gun
2. Click AddButton to add a new projectile
3. Select System: AlH3
4. Set Mole Fraction: 2, Frequency: 12000, Start Step: 1, Stop Step: 12500000
5. Set Fractional coords box: to 0, 1, 0, 1, 45, and 49
6. Enable Rotate
7. Velocity dummies: click AddButton, and select the top dummy atom followed by the lower one
8. Velocity direction: click AddButton
9. Set Energy: 0.05 eV, Atom Temperature: 450 K, Contact distance: 3 Å

To add the O2 projectile

10. Click the topmost AddButton to add a new projectile
11. Select System: O2 and repeat the above steps 4.-9. (reuse the previously added dummy atoms)
12. Change to Mole Fraction: 4 for the O2 projectile.


Except for Mole Fraction, the settings have to be identical for both projectiles as shown in the figure below



The Stop Step setting is used to realize an annealing phase at the end of the simulation

As a last step, we enable the molecule sink to remove any stray particles from the model.

1. Switch to Models Molecule Sink
2. Set Formula: * to remove all stray particles
3. Set Fractional coords box: to 0, 1, 0, 1, -150, and 150
4. Set Frequency: 12000 and Start step: 1

We are now ready to start the simulation using File Run. For the readers’ convenience, we provide the input file for this calculation as a download Ge200O25_O2_AlH3_fbMCCVD.ams


Due to the nature of the simulation, the calculation will run for an extended time. Depending on the hardware available, several days of runtime should be anticipated.

Results: Alumina Film on Ge(001)

After the simulation, we obtain the model structure shown in the figure below. The original surface is covered by an amorphous interface layer followed by a thick but loose and alumina film with many absorbed water molecules and OH groups as well as several chemisorbed aluminum complex moieties on the surface of the film.


We conclude this tutorial by noting several observations from the results:

Deposition Process
During the CVD simulation, the projectiles can be observed to appear just above the surface, getting absorbed to either react or desorb later, or getting reflected immediately. On average, more projectiles end up bound to the surface, which creates an amorphous alumina film growing over time.
Energy and Temperature
Consequently, the total energy of the system is continuously lowered as more chemical bonds are formed. In the beginning of the simulation, we observe strong temperature fluctuations, which gradually decrease over time. This latter observation is the result of an initially small model whose temperature can be easily affected by impinging particles.
Model Size
As this tutorial includes lengthy calculations, we opted for reducing the model size. In scientific studies, the use of a deeper and wider slab model might be more appropriate.
Alumina Film
Depending on the exact simulation condition (temperature, particle flux, particle kinetic energy etc.) films with different properties emerge. E.g. increasing the temperature generally reduces the water content in the film while increasing its density. The mole fraction of the flux of particles can affect the deposited film too. To obtain a clean alumina film, it is advisable to use an excess of oxygen during the deposition.
Ge(001) surface and interface
We generally observe the germanium surface to remain intact without any of its atoms desorbing. The interface between the germanium and alumina layers forms from an GeAl alloy layer, which in turn originates from the interaction between the AlH3 particles and the clean surface. With increasing oxidation, the deposited film separates into a mostly pure alumina film and an GeAl allow interface layer. Depending on the exact conditions, small alloys clusters are also be observed to move on top of the film during the phase separation.