TEMPERATURE-DEPENDENT PROPERTIES, DIFFUSION AND GROWTH ON THE SI(100) SURFACE, USING THE CAR-PARINELLO METHOD
Jerry Bernholc, NCSU
Scientific Significance:
The Car-Parinello method which combines molecular dynamics and density functional theory can be used to optimize geometries, find saddle points, and perform relatively short ab initio molecular dynamics simulations. Our current project involves ab initio simulation of high temperature motion of a Si adatom across a step edge on a Si(100) surface. The motion of such adatoms should result in the "step-flow" mode of growth, which is the most common growth mechanism in Molecular Beam Epitaxy. The overall goal of our project is to understand, in atomic detail, the growth modes of semiconductors and to identify the conditions for high quality layer by layer growth.
The Si(100) surface is technologically the most important of the Si surfaces, since this is the preferred substrate for the manufacturing of devices. The simulations utilized a stepped geometry, using a supercell that contains seven layers of Si atoms and one layer of H atoms, for a total of 232 Si and 72 H atoms. For studies of diffusion, three adatom have been deposited on this surface. We have performed a set of heating up "experiments" for this supercell. It is well-known that the Si(100) reconstructs to form a 2x1 dimer pattern, and that this pattern is stable at high temperatures. At low temperatures, STM imaging shows that the dimers buckle, forming a 2x2 or a 4x2 structure. Our simulations reproduce this structure well, and also show that the preferred reconstruction pattern depends on the kind and number of steps present at the surface.
It has been postulated that the dimers become effectively symmetric and/or disordered at higher temperatures, due to flipping of the dimer structures. However, our simulations, which were carried up to very high temperatures, show that the buckling of the dimers is still preserved, and that dimer motion is correlated over several dimer lengths at least. The details of the runs also provide information about the distribution of phonon frequencies and the softening of the surface phonon modes with temperature.
The current simulations focus on adatom diffusion on the stepped surface. This information is needed to simulate, via Monte Carlo methods, high temperature growth of Si structures. In this regime Si grows via the so-called "step-flow" mechanism, in which adatoms that have condensed on the flat regions of the surface (terraces) diffuse to the step edges, where they are more readily incorporated into the crystal. This mechanism results in the desired layer-by-layer epitaxial growth. Diffusion is a relatively slow process, since the adatoms need to cross a potential barrier and the "valleys" for adatom diffusion on Si(100) are particularly narrow, due to significant reconstruction.
The relevant time-scale is thus several picoseconds. The simulations so far have progressed to slightly over one picosecond. Although several aborted crossings were observed, no adatom jumps have yet occurred. We are planning to continue the runs to at least 10 ps, which should be sufficient to see several diffusion events. The results will serve as inputs to a Monte Carlo simulation of growth, which will be carried out over a much longer time period.
Numerical Approach and Performance:
The method combines an iterative solution of the electronic Hamiltonian with molecular dynamics. Since it is iterative, there are relatively large variations in convergence properties, depending on the complexity of the problem. For a perfect or almost perfect crystal, the electronic coordinates converge in 25-50 time steps, while a simultaneous optimization of electronic coordinates and atomic positions may require a thousand or more steps. (This is particularly true for highly distorted geometries and for metallic and nearly metallic systems.) In molecular dynamics or simulated annealing, the number of time steps is of course directly proportional to the length of the run. The time steps in the true MD simulation must be significantly shorter, compared to pure geometry-optimization, in order to preserve the constants of motion over a large number of time steps. Typical time steps in this approach in fairly complicated geometries range from 5 to 40 atomic units.
This code has a performance of 693 Mflops on a single processor of the C90, and 2.3 Gflop/s with 4 C90 CPUs. When the number of CPU is increased, the performance saturates rapidly with large overheads. The code has been ported to Cray-T3D. It runs heterogeneously on T3D with a C90 CPU as a front end. The electronic wavefuction parts are completely ported to T3D, which accounts for more than 95% of the running time of the code.
Back to Contents Page