Using coordination number to monitor a molecular dynamics simulation of  molecular beam epitaxy in a Lennard-Jones metal

Mathieu Bouville

 

 
 
 
 
 

A simple model of  molecular beam epitaxy (MBE)  of a Lennard-Jones metal [1] had been  simulated with molecular dynamics. The coordination number of the adatoms is estimated from the potential energy and used to determine if a given atom is in the bulk, in the top layer, on the surface as a monomer, a dimer, a trimer, etc.  The value calculated will depend on the cut-off radius (rc). Simulations with different cut-off radii in the simulation and during the calculation of  the coordination number are compared to determine the sensitivity of this technique to the choice of cut-off radius. 

 
 
Chime Plugin needed to view the below visualization

Movements of adatoms at the surface of a FCC crystal

Introduction

The aim of this project is to simulate the growth of a metallic crystal by molecular beam epitaxy [2] using the Lennard-Jones (LJ) potential .

. Atoms of the same species as those of the substrate travel toward a perfectly crystalline substrate with low energy.  They reach the surface but do not penetrate as in an ion beam. This study will serve as a stepping stone toward the study of MBE of semiconductors.  The simulation aims at implementing the MBE process before implementing a covalent potential (Stillinger Weber three-body potential [3] or Tersoff bond-order potential [4, 5]). 

Computational procedures

The substrate is FCC.  It is 5 unit cells long and wide and 4 unit cells thick.  There are 400 atoms in the substrate. Masses, and the coefficients s and e in LJ potential are all equal to 1 in LJ units. Periodic boundary conditions are used in the directions perpendicular to the growth direction to simulate a large substrate. Along the growth direction the top surface is free and the atoms of the bottom layer cannot move vertically in order to simulate a substrate whose bulk is infinitely stiff. The temperature is controlled through a Nose-Hoover thermostat. The pressure is held at 0 using a Parinello-Rahman barostat. The MBE process itself consists of three atoms impinging on the substrate with a low randomly directed velocity.

Quenching and data-base building

The coordination number of each atom or a “pseudo-coordination number” will be extracted from the simulations. As the coordination number only exists in perfect crystal we need another concept for this far-from-equilibrium process. The “pseudo-coordination number” Z’ of an atom will be calculated from its potential energy normalized by the one in the bulk: where Z is the coordination number in the perfect bulk crystal. Z’ will be called simply “coordination number” in the remaining of the text. In the bulk Z’ = Z (=12 in the case of FCC), in the top layer Z’ is expected to be a bit more than Z/2, for an atom on the top layer it will be less than Z/2 and it will be 0 for atoms which have not reached the surface yet. Hence it is possible to know if an adatom reached the surface, if it is in a dimer, a trimer, if it was buried by later-arrived adatoms, etc. As the potential energy depends on the cut-off radius, so will the coordination number. 

In order to use the coordination number to monitor the MBE process typical coordination numbers had first to be known. Of interest were bulk atoms, surface atoms, monomers, dimers, aligned and triangle trimers. Quenching simulations allowed to calculate coordination numbers at 0K: atoms are put on the surface of a substrate in dimers, trimers, etc. and the temperature is decreased to 0. Results are given in the following table.
Cut-off =1.6
Cut-off =2
Cut-off =2.5
potential energy
coordination number
potential energy
coordination number
potential energy
coordination number
Bulk
-9.37
12
-12.72
12
-14.8
12
Surface
-6.25
8
-8.6
8.1
-10.95
8.85
Monomer
-3.1
3.98
-4.15
3.9
-5
4.05
Dimer
-3.88
4.98
-5
4.8
-6
4.85
Aligned trimer
-4.66
5.97
-6
5.65
-6.98
5.65
-3.95
5.05
-5.05
4.77
-6
4.87
Triangle trimer
-4.66
5.97
-6
5.65
-6.97
5.65
-3.95
5.05
-5.3
5
-6.21
5.03
Potential energy and coordination number of different types of atoms for cut-off of 1.6, 2 and 2.5

Data Analysis

Simulating at different cut-off radii (1.6, 2, 2.5) and analyzing with the same cut-off as in simulation turned out to be useless as the trajectories were very different (for instance there was no trimer for rc = 2). 

Therefore one simulation was run (with a cut-off radius of 2.5) and the coordination number was then calculated from this unique simulation with three different cut-off radii (rc = 1.6, 2, 2.5). The figure shows the coordination numbers of the three adatoms for rc = 2.5. The results look the the same with the other two radii). We can see successively monomers 401-402 dimer, trimer centered on 401 and 400-401 dimer. The thermal noise as expected does not allow to know if the trimer is an aligned or a triangle one. Quenching may solve this problem but in that case the coordination number cannot be used to monitor the simulation while it is running.


Graph of mean square displacement vs. time taken from simulation data

Conclusion

Coordination was used to monitor molecular dynamics simulation of the displacements of adatoms on the surface of a perfect crystal as in MBE process. It was found useful to recognize bulk atoms, surface atoms, monomers, dimers, and trimers. However more subtle differences (between aligned and triangle trimer) are hidden by the thermal oscillations of the atoms. The cut-off radius, in the analysis, was found to have little qualitative importance: short and long cut-off enable to recognize the basic patterns but at non-zero temperature none is able to see more subtle differences.

References

[1] J. E. Lennard-Jones, Proc. R. Soc. London A 106, 463 (1924).

[2] see for example E. H. C. Parker (editor), The Technology and physics of molecular beam epitaxy (Plenum Press, 1985) or R. F. C. Farrow (editor), Molecular beam epitaxy : applications to key materials (Noyes Publications, 1995). 

[3] T. A. Stillinger and F. H. Weber, Phys. Rev. B 31, 5262 (1985). 

[4] J. Tersoff, Phys. Rev. B 37, 6991 (1988). 

[5] J. Tersoff, Phys. Rev. B 38, 9902 (1988).