Calculation of the Electronic and Optical
Properties of Nanoscale Systems.
A thesis submitted to Dublin Institute of Technology
For the award of PhD
Dr. J. D. O'Mahony, School of Physics and FOCAS, DIT
Dr. J. P. Lewis, Department of Physics, WVU.
At the nanometer length scale, the size of surface features in crystalline semiconductor systems is of the same order as the electron wavelength. This can result in unusual behaviour in the systems electronic, magnetic and optical properties due to electron confinement effects. Such effects can have practical and commercial applications and are currently the subject of considerable study in the disciplines of theoretical, computational and materials technology within nanoscience. This thesis uses molecular dynamics computational methods to examine such effects in the electronic structure of semiconductor-based crystalline systems. Three unique surfaces were studied in detail - the SiC(111) surface, the SiC(100) surface, and the prototypical In-Si(111) surface. Silicon carbide is of importance in the development of semiconductor technologies due to its physical robustness and relatively high power capabilities. An understanding of surface metallisation in semiconductors is of paramount importance since modern technology relies on the interaction of metals with semiconductors in integrated circuit and device construction. If Moore’s Law is to be adhered to, transistors must become smaller and the metal contacts between transistors must likewise shrink. This work explores the possibility that potassium deposited on the SiC(100) surface may provide a solution for nanoscale contacts between devices on this surface. Using modified and highly efficient molecular dynamics code, the energies and reconstructions of a number of possible surface configurations were studied in detail, resulting in proposed new candidates for surface reconstruction for a range of coverages of potassium on the SiC(100) surface. The SiC(111) surface has previously been shown to undergo an interesting metal-insulator transition where the surface band states split. This has been observed by experimentally probing the surface states with scanning tunneling spectroscopy and photoemission techniques. By applying ab-initio molecular dynamics techniques to simulate this surface, this research has found compelling evidence for the actual mechanism that results in this transition. A number of time-dependent simulations of the surface in question were carried out, over ranges of tens of thousands of picoseconds. The results show that the surface is dynamical in nature. Furthermore, the transition is shown to be due to a soft phonon interaction on the surface, and thus surface dangling bonds are seen to split because they are in constant motion. Finally, computational studies of the In-Si(111) surface are also presented. The results indicate a dynamical surface exhibiting surface phonon effects, similar to the SiC(111) surface studied and metallisation in a similar vein to results obtained for the K-SiC(100) surface. The study of the In-Si(111) surface therefore represents a natural progression in studies of this nature. The computational work presented here was carried out using the FIREBALL suite of tools. During the course of this study, the codebase was rewritten and modernised to improve performance and to allow for easier future modification. The extensive changes to the code are discussed, as are its potential future applications in the field of computational solid state physics. Practical methods are presented that allow for the work to progress to the calculation of optical transitions directly in FIREBALL, with a full description of how a reflectance anisotropy spectrum could be calculated as a logical extension of the present work. The calculation of a reflectance anisotropy spectrum would be of considerable interest to experiments in the field.
In Solid State Physics unusual behaviour can be seen in a material’s properties when nanometer dimension scales are realised. Novel behaviour can manifest itself in the electronic or optical properties of the system and may lead to technological breakthroughs when these new properties are exploited. An example of such applications are quantum well lasers which use electron confinement to achieve shorter wavelength emission and are more efficient than conventional laser diodes. At short length scales interesting emergent behaviour such as quantum wire characteristics, Mott-Hubbard transitions or unexpected surface metallic effects are observed. Understanding these properties is crucial to further development of physics. As technology approaches smaller and smaller designs, an understanding of the novel physics at this level is required so that such behaviour can be exploited for new applications.
When attempting to control surface properties to engineer nanoscale structures for technological applications, the experimental costs can be very high. For this reason much work has gone into the development of computational tools over the last number of decades. In 1965 theoretical physicists, Walter Kohn and L.J. Sham published a paper on a then new technique called density functional theory (DFT). At the time they felt that this would not have a large impact on the academic community. However it opened the door to pure ab-initio molecular dynamic and electronic calculations for systems of complexity beyond the hydrogenic system. Until then, solving anything more complicated than the hydrogenic system by analytical or numerical methods was an arduous task, and thus mathematical modeling of nanoscale systems was not viable.
With the advent of modern computers and developments in computational methods, it is now commonplace to simulate systems of hundreds or thousands of atoms. Many different possible permutations of atomic species as dopants, or as surface adatoms, can be explored before physical experiments are carried out. This allows for more efficient use of experimenter’s time and increases the predictability of such experiments. Using computational tools, a comparison can be made between theoretical and experimental results or the predictability of experiments can be evaluated before actually carrying out more expensive work.
Computational tools also allow a researcher to determine the mechanism by which some previously unexplained experimental results have arisen. For example, in the Sn/Ge(111)-√3 × √3 to 3 × 3 reversible phase transition (RPT), computational work has been used to reliably describe the mechanism responsible. Experimentally, the mechanism was contentious for some time before a computational study of the problem showed that the interaction of soft phonons, which could not be resolved in the time-domain by current technology, were present and are most likely the cause of the RPT.
Another advantage of computational over experimental methods is the cost. All fundamental research is expensive, but in using computational methods the expenses are in computers and programming. With experimental work the production of raw material of sufficiently high quality and the generation of extreme isolated environments, for example ultra-high vacuum systems, is a substantial investment. Computational methods allow for a large number of “experiments” to be carried out before such investment is necessary for real-world experiments. This work discusses the development of a computational physics package as well as work carried out, using computational methods, to explore some interesting behaviour of nanoscale systems.
In Chapter 2, a number of techniques in computational solid state physics are discussed. Some approximations can only be applied to very small systems, in terms of number of atoms. This also is the case with analytical quantum mechanics methods. Other approximations lose accuracy with respect to expected experimental results, but allow for simulations of far larger supercells. A supercell in computational solid state calculations is an ideal repeating unit of, for example, some bulk crystal.
There are a number of methods not covered by this thesis, such as quantum Monte Carlo (QMC), where random numbers are used to evaluate the outcome of quantum mechanical operators or semi-empirical methods. Semi-empirical and empirical methods take some values for parameter fitting to run simulations, as a result of this they cannot evaluate electronic properties of a system.
In Chapter 3, the FIREBALL package used in this work is extensively discussed with particular reference to the uniqueness of the formulation used and the linear combination of atomic orbitals (LCAO) method. This method allows for the precompilation of tables of interactions before performing any simulations and allows FIREBALL to simulate far larger systems than other ab-initio MD packages for the same amount of computational usage. The significant enhancements to FIREBALL made during this project are also discussed.
Chapter 4 discusses the generation of FIREBALL’s precompiled tables- the ‘Fdata’, and how this is carried out within the current FIREBALL package. Chapter 5 discusses in detail the rewriting of this program to create the new package FIREBALL-Lightning, which is currently a non-MD package. The FIREBALL-Lightning package will ultimately be further developed to succeed FIREBALL.
Chapters 6, 7 and 8 discuss the actual experimental work carried out for this project. This work was carried out using both FIREBALL and FIREBALL-Lightning, with both a new and an older implementation of the Fdata generation program. The metallisation of the silicon carbide 3C (100) surface was explored in order to better understand the unique properties that were first shown by experimental methods.4 Silicon carbide receives much interest in the literature as it may be a viable substrate for semiconducting technologies in harsh environments due to its robustness and high power stability. The metallisation of this surface shows that even at low coverages the surface may be able to conduct charge, which is of paramount importance in device design.
Chapter 7 deals with another surface of silicon carbide, namely the (111) surface. The “silicon-not-so-rich” surface has been shown experimentally to exhibit a ‘Mott-Hubbard transition’, which is a phenomenon whereby a surface band is seen experimentally to split into two bands, altering the conducting nature of the surface. This is studied in detail and an alternative explanation for the system’s behaviour is proposed. A greater understanding of this phenomenon is required as length scales shorten in industrial fabrication and design.
Chapter 8 deals with the so-called “quantum-wire” behaviour of the In-Si(111) surface. Current experimental and computational results are discussed in depth and preliminary work towards developing a system to calculate a dynamic reflectance anisotropy spectrum (RAS) from first principles is discussed. This surface presents a reversible phase transition with temperature which may be due to a soft-phonon effect. This would mean that, in order to make a genuine comparison of computational optical spectra (such as RAS), not only is a ground state evaluation of the system required, but also a dynamic range of further states would be needed, due to this soft-phonon effect. This may well prove the key to further understanding this controversial surface.
There are many different levels of theory for simulating materials. Generally the more exact the method, the more computationally expensive it is. At the simplest level there are empirical methods. Empirical methods are fitted to experimental results via various parameters. New knowledge is then gleaned about the experimental results from fitted simulations. Such methods are very popular in the modeling of highly complex systems that consist of many atoms, such as biological systems. In most empirical methodologies the equations of motion are invoked so the Hamiltonian is constructed as in classical mechanics, such as:
and forces acting on particles within the system by:
and the Hamiltonian for molecular dynamics would consist of something like:
where the above symbols have their usual meanings and the constant, k, is the parameter matched to experimental results. In this example, the new information gained might be the time-development of position by using the equations of motion. Because of the low computational cost, empirical methods are used mainly in computational biology and in modeling of weakly-interacting systems. Such systems consist of thousands of atoms and as a result are outside the limits of current computing power for other methods. The advantage of empirical methods is that they are inexpensive and, once parameters are tweaked, numerical evaluation is relatively simple. The disadvantage is that we have no knowledge about the electronic structure of the system. What is gained is some insight into the time evolution of the system, i.e. the molecular dynamics (MD).
One of the earliest references to molecular dynamics was originally outlined by Adler and Wainwright in 1957 using classical models, namely Newtonian forces, to calculate the positional movements of solid spheres, which represented atoms in a liquid. The principle is that computational models are used to simulate interactions between all atoms, as if they are moving continuously. MD allows for a study of what is going on within a system (such as a semi-infinite crystal, a nanocluster, a biological molecule or system derived from these components) at specific temperatures or under other conditions such as varying pressures, etc.
More exact methods use a quantum mechanical approach to calculate the forces. These are known as ab-initio methods. Ab-initio, from the Latin for “from the beginning”, conversely, has no fit parameters. The fundamental building block of all ab-initio methods is the solution of the Schrödinger equation:
where H is the Hamiltonian and ∕epsilon represents the energy Eigenvalues for the Eigenvectors Ψ. This implies that computational simulations can predict the outcome of experiments that have not yet been carried out. It also provides an effective way to evaluate the electronic structure of the system. Forces, electron behaviour, and a wealth of information can be determined from such simulations. Because ab-initio calculations are completely independent of experimental physics results, comparison of the results from the two separate approaches can give further insight into the mechanisms at play in any system of interest. The disadvantage of ab-initio methods is that, comparatively, they are computationally far more expensive. Some implementations cannot simulate more than a few hundred atoms or less, due to the limitations of current computing technologies. This will be explained in greater detail in section (2.2).
Density functional theory (DFT), originally formulated by Thomas and Fermi in 1927, but taking its modern form under the work of Kohn and Sham in 1965, is a method for solving the Schrödinger equation for use in MD and electronic structure calculations. It is based on the theory that the many electron problem can be substituted by an electron density projected over space. This approach therefore reduces the number of dependent variables in the Schrödinger equation to just one. The Kohn-Sham approach, and the Harris-Foulkes method utilised in FIREBALL, requires a number of approximations to calculate the forces involved. These approximations and DFT are discussed in section (2.2.1) to (2.2.4).
Tight binding theory, which is similar to the linear combination of atomic orbitals (LCAO), was originally introduced in 1929. The implementation of this in FIREBALL as well as in other popular MD packages is discussed in this chapter and the next. As a basic summary, the electron wavefunctions are modeled as being centered on the atom to which they belong, and allowed to go to zero at some cutoff distance from the atom. This approximation is acceptable as it can be shown that the energy Eigenfunctions of any bound electron fall away quickly with distance from its host atom. A more in-depth discussion of this methodology is presented in the rest of this chapter.
In Chapter 5 an outline of the implementation of exact exchange within the FIREBALL suite is presented. As discussed in Sections (2.2.1) and (2.2.2), exchange and correlation are fundamental issues in computational physics. Exact exchange is a method for modeling the exchange in a more accurate way than was previously possible. This ongoing work is developing towards a platform for calculating a reflectance anisotropic (RA) spectrum from first principles which it is intended to ultimately merge with the current rewriting of the FIREBALL code using the more modern Fortran 90 language. Fortran 90 allows for further optimisation of the algorithms used. These are also discussed in Chapter 5.
Starting from the time-independent Schrödinger equation for an isolated N-electron system, and using the Born-Oppenheimer non-relativistic approximation given by:
where Ĥ is the Hamiltonian operator,
in which T is the kinetic energy operator, VNE is the electron-nucleus attraction energy operator and VEE is the electron-electron repulsion operator. In more explicit form, each of these terms can be rewritten as:
where all symbols have their usual meanings and v is the potential acting on electron i, due to nuclei of charges Zα, given by:
In equation 2.2, E is the electronic energy and Ψ = Ψ(x1, x2, x3, ... xn) is the wavefunction and xn represents the coordinates of each, sometimes with the addition of the spin as a fourth dimension.
Equation 2.2must be solved within appropriate boundary conditions: the wavefunctions must be well behaved throughout all space, and decay to zero at infinity, or be subject to some other periodic boundary condition for a regular infinite solid.
The Born-Oppenheimer approach separates the Schrödinger equation into an electronic part and a nuclear part. It is the first approximation made in almost all quantum mechanical calculations, as it allows the wavefunctions to be separated into the kinetic energy, the electrostatic potential due to the nucleus and the potential due to other electrons. Of course, in the presence of other fields, there would be more terms in, for example, equation 2.2. The assumption is that the electrons respond immediately to nuclear motion, and this can be dealt with separately. With systems of more than one electron, it greatly simplifies finding solutions to the Schrödinger equation.[15, 16]
Early attempts to numerically solve the Schrödinger equation for many-body systems were based on Hartree’s method. Hartree was reportedly seeking a non-empirical method to solve the Schrödinger equation for complex systems by assuming that any one electron moves in a potential that is the spherical average of all other electrons and the nucleus. The spherical average is assumed to be centered on the nucleus and calculation of it by numerical integration for any one electron requires knowledge of the wavefunctions of all other electrons and the nucleus. The method requires an approximation to be made which is then used to solve for the wavefunctions of all bodies, whereupon an improved approximation is arrived at, and the process repeated from this new starting point. The method became known as the self consistent field (SCF) theory because of this process of iteration. Later, the method was clarified by Fock with an explanation of the validity of its operation using the variational principle. The variational principle essentially states that any Ansätze to the ground-state solution of a system must have higher energy than the exact solution. This leads to an iterative approach to solving the Hartree-Fock equation numerically.
The variational principle states that a measurement of some observable, for example energy, is always higher than or equal to the actual ground state energy of a system (that is, it is never lower). Given that:
where the square brackets indicate that E is a functional of ψ, and all other symbols are as before, any measurement of E results in a value that is, by definition, greater than or equal to E0. Therefore:
Thus the energy computed from a guessed Ψ is an upper boundary to the ground state energy. The true ground energy can be calculated by finding the minimum of the functional E[Ψ]. This is known as the variational theorem, which is very useful in these calculations.
Using the variational theorem, we can now look for a method to solve for the ground state of a system by applying an approximation to Ψ and E0 and iteratively converge on a better solution.
The Hartree-Fock (HF) approximation is a method whereby the wavefunction, Ψ, is approximated as an antisymmetrised product of orthonormal spin orbitals Ψi(x). Each of these spin orbitals is itself a product of spatial orbitals ∕phik(r) and a spin function σ(s) = α(s) or β(s). The method essentially is the determination of these orthonormal orbitals in such a way as to minimise the energy as calculated in equation (2.2.1).
So the Hamiltonian for any system is the sum of smaller Hamiltonians that make up the system:
By normalising the integral of the HF wavefunction, as done in all real-world quantum mechanical calculations, we can develop a method to calculate the expectation value of the energy. The expectation value of the energy is then given by:
These integrals are all real, J and K are positive. J is known as the Coulomb integral and is the Coulombic repulsion between electrons. K is the exchange integral, for which there is no classic counterpart. It comes about because the electrons are two identical particles and the integral can be seen as the interaction between two states in which the coordinates of the particles are exchanged. The exchange integral is therefore a result of the Pauli Principle. For more detail on this see ref [ sec. 1.3]. The solution must proceed iteratively, since the orbitals that solve the problem also appear in the operator. Hence the name, the set of equations become self-consistent. This method gives orbitals that are highly localised, and are useful for some purposes, such as when evaluating a highly ionic system.
The usual implementation of the HF theory, and other related methods, is to employ some set of one-electron basis functions, for which orbitals are expanded and many-electron wavefunctions are expressed. This allows the problem to be transformed into one or more matrix Eigenvalue problems of high dimension, where matrix elements are calculated from arrays of integrals evaluated for the basis functions,χp, for example, see equation 6 of ref [ sec. 1.3].
Up to this point, a general framework for solving the functionals involved in these sort of calculations has been presented. However, this is limited in a number of ways. For example, any purely quantum mechanical effects have not yet been accounted for, and certain modifications now need to be made.
One of the inherent errors of the HF formulation is that it does not take into account the so-called “correlation energy”. To be precise, the difference between the exact energy of a known system and the energy calculated by the HF method is known as the correlation energy, i.e:
This difference comes about because within the HF, the wavefunction is approximated by a single Slater Determinant, or some other approximation. This causes an error because the exact solution is never that simple.
There are many approximations in the literature to deal with the correlation energy and to correct for this error. Often it is bundled in with the exchange energy and both exchange and correlation are dealt with simultaneously, called the exchange-correlation (XC) energy. Popular methods for dealing with the XC are discussed below in implementation of Density Functional Theory.
The final fundamental piece of theoretical information required for understanding the basic background theory is the use of electron density. In terms of Ψ, the electron density ρ(r), which is the number of electrons per unit space of a given state, is
By inspection, it can be seen that this is a non-negative function and that the integration over all space will be equal to N, the total number of electrons in the system.
The electron density is very useful for numerical calculations, and can be combined with the HF or a similar theory. It allows for the rewriting of the Schrödinger equation in terms of a single variable, namely ρ(r), which is far easier to work with. For the many-electron Schrödinger equation, the T and Vee terms are dependent on the overall number of electrons in the system ( = ∫ ρ (r) dr ) and the Vne term is dependent on electron density in the vicinity of r.
An expression in ρ(r) is acquired by applying the variational principle to the electron density formulation of the Hartree Fock method (see Yang and Parr’s book on the subject, ref. []). The origin of the exchange-correlation error discussed above is in the Vne term and approximations to solving this issue are explained in Section 2.2.5.
For an N electron system, the external potential, v(r) completely fixes the Hamiltonian when considering the ground state only. Thus, N and v(r) determine the properties of the ground state of the system. Hohenberg and Kohn proved in 1965 that “The external potential v(r) is determined, within a trivial additive constant, by the electron density ρ(r).”
Since ρ(r) determines also the number of electrons, it can be shown that ρ(r) also determines the ground state of the wavefunction and all other electronic properties of the system. Consider the electron density ρ(r) for the non-degenerate ground state of some N-electron system. N is determined by simply integrating over all space. However, ρ(r) can also be used to determine v(r) and, therefore, all other properties of the system.
The issue of exchange correlation then arises because the external potential term will always result in a resultant E that is higher than the actual E.
The Hamiltonian can be rewritten in terms of electron density, giving:
The exchange correlation energy can then be expressed as part of the Vee term as:
where J[ρ] is the classical electron-electron repulsion discussed above in section 2.2.1 and Vxc [ρ] is the nonclassical electron repulsion, or the exchange-correlation (XC) term.
Density functional theory takes the principles outlined above to a more workable form for numerical calculation. The theories of Hohnenberg, Kohn and Sham[8, 20] must be outlined here in order to clarify the solutions used.
Kohn and Sham succeeded in showing that if we take the energy due to the XC term as given by:
The corresponding potential is then:
Applying this to an atom, a molecule or a solid is essentially the same as if the XC energy for a non-uniform system can be found by applying the uniform electron gas results to infinitesimal portions of the non-uniform electron distribution, each having ρ(r)dr electrons, and then summing over all space. This is known as the local density approximation (LDA):
The LDA is applicable to systems with slowly varying densities but cannot be justified for highly inhomogeneous systems. It is most valid when applied to calculations on solids. The Perdew and Zunger implementation of LDA is used in the FIREBALL suite. The other XC approximation available within the FIREBALL suite is a Generalised Gradient Approximation (GGA), and the implementation of Lee, Yang and Parr ref [] is used in FIREBALL. The GGA takes what we have for the LDA but adds a gradient in the n terms, resulting in a better approximation to the exact result but at a computational cost. Rather than hypothesizing uniform density “pockets” of infinitesimal size, these pockets now contain a gradient. An example of such an implementation can be found in ref [].
Combining all the above into a formalism gives nearly equivalent results to solving for the system analytically, especially if the GGA is chosen over the LDA. However, the computational cost is still very high. The method listed above requires the modeling of every electron in the system, even those that have little or no effect on the chemistry, the core electrons. For this reason, we introduce the idea of pseudopotentials, as now explained.
Around the nucleus, the wavefunctions of electrons are not very smooth. In this “core” region, the electrons do not interact in any major way with the solid. When attempting to model the electron-electron interactions, a popular approximation is the use of “pseudopotentials” in this core-electron region. Essentially, the electron wavefunctions at this level are replaced with smooth pseudopotential wavefunctions. This is not an uncommon approximation in molecular dynamics implementations. The fact that the valence electrons are responsible for chemical bonding, whereas core electrons are generally unchanged in a different environment, makes this a very reasonable approximation. This means that only the valence electrons need to be accounted for in a self-consistent calculation. Computationally, the pseudopotential wavefunctions are far easier to deal with when calculating the electron density as outlined above. By using the pseudopotential approximation, the Hartree potential, vH, and the exchange correlation, vxc, terms are now only evaluated for the valence electrons in veff and only for the valence density, ρV. The core electrons are accounted for by replacing the external potential, vext with a pseudopotential, vextPP . Hence:
There are a number of methods in the literature for modeling the core electrons as a pseudopotential. The main criterion is that at some point, such as the boundary between the core and valence electrons, the effective potential is equivalent for a single atom of the element in question. [24, 26]
There are many ways to calculate the pseudopotentials, and no one method is significantly superior in any way to another. Commonly followed prescriptions include Kleinman and Bylander, Vanderbilt, Troullier and Martins or Hamman and Schluter. The pseudopotential approximation is popular because it considerably reduces the number of electrons treated explicitly in the calculations, greatly decreasing computational cost. Another advantage is that the wavefunctions are smoother farther from the core, as the nucleus affects the core electrons far more than outer electrons. The model used to describe the remaining electrons is simplified by ignoring the effect of the core.
There are many differing MD and DFT packages that have emerged since the first formulation of the approach. The differences between them are in some cases quite slight. One major implementation is the GTO, the “Gaussian Type Orbitals”, or Gaussians. This models the integrations of Gaussian functions which can be reduced by way of the Gaussian product theorem to fewer and fewer required integrals. Gaussian functions are used to model the orbitals in an LCAO implementation, where a four-center integration can be reduced to two centers and then to one. It is based on the work of S.F. Boys in 1950.
Another major implementation is the plane-wave method (PW), as used in the Vienna Ab-inito Simulation Package (VASP), ABINIT and CASTEP. Within PW DFT, it is usual to employ a pseudopotential scheme similar to that previously discussed.[33, 34, 35] The essential difference is that the wavefunction of the entire system is modeled by a “plane-wave basis set” in which, unlike in an LCAO implementation, the wavefunction is constructed to model the entire system. The basis functions are then orthogonal, which is lacking in the FIREBALL implementation until the application of the Lowdin transformation. The PW method makes use of a Fourier transform of the overall electron wavefunctions to model the electronic behaviour of a system. The limitation of this method is that it becomes very expensive computationally with increasing cell size, whereas FIREBALL becomes more expensive with the number of species involved.
In this chapter the background theory to DFT is discussed as well as popular approximations and implementations. Like FIREBALL, the SIESTA package uses a tight-binding model based on the Sankey-Nikwelski method. This is discussed in greater depth in the next section. Unlike PW, the tight-binding model is based on spherical orbitals which are atomically centered. This is explained in depth in the next chapter.
This chapter deals with the many aspects of ab-initio calculations specific to the FIREBALL suite. There are a number of programs that make up the entire package, and these are dealt with as required. In Chapters 4 and 5, the Create and Create-Lightning programs are discussed. The Create program generates a number of tables which contain the data used by FIREBALL to carry out its simulations. These tables contain all integrals computed on a numerical grid for all one-, two- and three- center interactions. These tables are collectively known as the Fdata. During atomistic simulations the value of the integral at a specific point is gleaned from the Fdata tables at runtime. The Create program uses pseudo-atomic wavefunctions for each the of species of interest which is generated by another program, Begin. The Begin program was not developed by the author and is only mentioned where required.
The FIREBALL package, which was used in all computational work in this project, uses an ab-initio tight-binding formulation that was first developed by Sankey and Niklewsky. It is based on norm-conserving pseudopotentials and DFT within the LDA or GGA. Rather than the Kohn-Sham functional, the Harris-Foulkes functional[13, 38] is employed. The main difference between the two is that the Harris-Foulkes functional depends solely on the input charge density, whereas Kohn-Sham is self-consistent. Self-consistency, as in the Hartree-Fock method, is iterative, being dependent on both input and output electron densities. FIREBALL also makes use of slightly excited pseudo-atomic orbitals within a localised basis set for the atoms to calculate the charge density. The electronic Eigenstates are expanded as a linear combination of these pseudo-atomic orbitals. They are slightly excited due to the boundary condition that they go to zero at and beyond some radius, rc, that is:
which confines the wavefunction, and thus raises the energy levels, “exciting” the electron. These orbitals are referred to as “fireballs”. When analytically solving the atomic problem, the wavefunctions extend to infinity.
In equation 3.2.1, the first term is the band structure energy, and is the sum over the occupied Eigenstates, ϵn, of the effective one-electron Hamiltonian:
The potential V,
is the sum of the ionic potential, vion(r), normally represented by a pseudopotential, a Hartree potential (that is, the classical electron potential in its ground state), and the exchange-correlation potential, Vxc. The rest of the Harris Foulkes equation is broken down as the average electron-electron energy, Eee:
the ion-ion interaction energy Eion-ion
where Z is the nuclear (or pseudopotential) charge of atom i at position Ri, and the exchange-correlation energy, Exc.
MD is performed by evaluating the forces, Fi, on every atom i where:
The previously-mentioned pseudoatomic “fireball” functions are used to solve the Schrödinger equation. However, the “atom in a box” boundary condition of the wavefunction raises the electronic energy levels as a result of the confinement. The cutoff radii in this case are carefully chosen to minimise errors with respect to the free atom state and so that the electronic Eigenvalues remain negative. The method used for choosing such cutoff values is discussed in depth in reference [].
Using these fireballs to evaluate the total energy, the forces on an atom at position Rl are determined by taking the derivative of the total energy with respect to Rl. The band structure force is evaluated using a variation of the Hellman-Feynman theorem.
Energies are found by solving equation 3.2.1. Using ρin as the input density, the sum of the combined spherical atomic-like densities of all electrons in the system in:
where ψi(r - Ri) are the FIREBALL wavefunctions used to solve the Schrödinger equation. The value of ni is the occupation number which is the number of electrons occupying that specific atomic-like spherical density. The total energy is then evaluated by using the a reference “atomic density” by taking ni to be the occupation of these pseudo-orbitals in the neutral atom rather than by self-consistent methods. This has been shown to yield results that are more accurate than the Kohn-Sham equations when not used self-consistently, and to be relatively inexpensive computationally.[37, 38]
The FIREBALL boundary condition offers two other benefits to the computational evaluation of molecular problems. Firstly, very sparse matrices can be generated for large systems because the range of hopping-matrix elements on different atoms is limited. This achieves a boost in computational efficiency. Secondly, somewhat serendipitously, the slight excitation of atoms yields a good representation of solid-state charge densities because it accounts for Fermi compression in solids.
The exchange-correlation options available in FIREBALL are uniquely suited to large-scale ab-initio MD. There are currently three implementations. However, one of them is largely superseded by another, and remains for legacy purposes. These three are the Sankey-Nikwelski method, the Horsfield method and the Generalised Sankey-Nikwelski method, which is also known as McWEDA or OLSXC.
The Sankey-Nikwelsky method only applies to a minimised basis set (specifically sp3). It has also been called “nearly uniform density approximation”. It is based on the idea that the matrix elements of the exchange-correlation potential can be rewritten as:
where the μ and ν indicate some quantum number arrangement and
is termed the “average density”. An approximation is required when the divisor, above, is equal to zero, and this was done for the sp3 basis set. This makes it unsuitable for any problem with a significant valance-band electron density, such as transition metals.
The Horsfield approximation is well suited to molecules and clusters of atoms. It is based on a many-center expansion of (ρ(r) = Σ(ρi (r)), above. By treating the problem as two separate cases, it is shown that:
known as the “atom” case in FIREBALL, and
known as the “ontop” case.
In each of these cases, the atomic sites of ψμ, ψν are given by i and j, respectively. In the “atom” case i = j, whereas in the “ontop” case, i ≠ j.
This has been shown to be an accurate approach in many cases. However in the atom case, there can be some discrepancy, and additional numerical integrations are required.[40, 43, 44] It is also quite expensive computationally.
The McWEDA approximation is well suited to problems based on solids. It combines both of the other two methods, as well as generalising the Sankey-Nikwelsky method.
To generalise the Sankey-Nikwelsky method, the atomic orbitals are replaced with spherically symmetric orbitals when calculating the ρ terms. This overcomes a number of issues. Firstly, we do not get a zero overlap for intra-molecule orbitals (e.g. two s orbitals on a single atom). Furthermore, this means that regions of positive overlap do not cancel out regions of negative overlap, so the “importance” of a sample is no longer skewed by such a measurement. This is done by separating the orbital into its radial and angular parts, taking the root-mean-square of the radial and replacing the angular part, i.e. the spherical harmonic, with that of a perfect sphere (i.e. Ylm goes to Y00 regardless of orbital type.
The next step in coming to the McWEDA method is simply to plug in the resultant expression, this time up to the second order, into the Horsfield method. This results in:
where, as before, μ and ν indicate some combination of applicable quantum numbers for the species in question and ϕ indicates the spherically averaged orbitals, and i and j again indicate the site of the orbital center and
for the “atom” case and
for the “ontop” cases.
The Harris-Foulkes functional yields excellent results for strongly covalent systems and has been shown to give results comparable to the LDA in non self-consistent calculations. However, it is limited in that unsatisfactory results are produced by the modeling of systems comprising species with large electronegativity differences. This is because in Harris-Foulkes we assume the occupation numbers of each orbital to be equivalent to the neutral-atom values and second-order errors are produced as an extension of this.[38, 41] In ionic, or even in non-pure-covalent systems, electrons move from their neutral atom configurations and the solution of these calculations requires an SCF approach.
The SCF reformulation of the Harris-Foulkes functional is referred to as DOGS, for the authors of the paper (Demkov, Ortega, Grumbach and Sankey in ref. []). Within the DOGS reformulation, a variance in the value of ni is introduced. Thus the energy functional is now dependent on n and we vary the input density by changing the n term to (n0 + δn) in eq. 8, where n0 is the neutral atom density, and apply a self-consistent calculation about this new formulation for population. This incorporates long-range effects, rather than the highly localised Harris-Foulkes non-SCF functional.
The final element of the FIREBALL method relies on the fact that the previously mentioned fireball orbitals are used. It has previously been shown1 that no more than three centers (orbital centers) are required to calculate all applicable integrals. As the FIREBALL method defines that any orbital goes to zero after the limit of distance of rc1 + rc2, (i.e. for a two-center interaction, given that all matrix elements are zero outside of the limits of their cutoff radii, the maximum distance that need be considered is rc1 + rc2), the calculations involved lend themselves to being pre-calculated and tabulated.
This very property enhances the speed of FIREBALL calculations greatly. It means that for all possible interactions of the species in any system of interest, based on the XC, on the cutoffs and on which species are involved or chosen, a table of integrals spaced on a numerical grid can be compiled before any molecular dynamics simulations are carried out. Within the MD simulation itself, the data in these tables is simply interpolated for the required situation from this spatial grid.
The benefit of this is that the integrals requiring tabulation can be pre-calculated en masse. The number of integrals required grows as order N3 with the number of species in the study. Thus pre-calculated integrals, which we call Fdata , lend themselves to this kind of processor-scaling and is especially beneficial in computational efficiency.
In the years since FIREBALL was first developed, the Fortran language in which FIREBALL is written has been updated and improved - FORTRAN77 became Fortran 90, which in turn has led to the minor upgrade that is Fortran 95 and the more heavily-revised Fortran 2003. One important new aspect introduced in Fortran 90 is variable control, which eliminates precision issues common in older versions of Fortran code, where double-precision may be declared for some variables but not others. Double-precision simply means usage of 16-bit real numbers instead single-precision 8-bit storage. Precision issues can also cause a problem in compile time. The new package removes these issues with an implementation that allows the user to decide on the requisite precision at compile time only.
Another Fortran improvement is the introduction of “derived types”, similar to “structs” in C++. This allows groups of variables to be sorted together in an encompassing variable folder (in the sense of a filesystem). Previously, there would be a large "library" of global variables to take care of information such as the atomic number of a species in any large program. Globally required variables can now be bundled together in derived types. For example, this means that under the type known as "species", a developer can access all the data required on that species, such as its number of shells, cutoff lengths, the name of the species and all other data pertinent to that data type.
A further development in the Fortran language allows for much better control of arrays. In older versions of the package, arrays had to be declared to a maximum possible size at compile time. Now they can be dynamically allocated, increasing the efficiency of the code. Arrays can also be operated on directly rather than iterating through the array one space at a time to carry out a summation or multiplication. This is a major simplification in writing code and allows for compiler- or architecture-specific optimisation of array handling.
The final notable newer feature of Fortran is the introduction of modular programming. A module is a self-contained unit of code for a specific purpose. The need to turn certain procedures "on" or "off" by input files has effectively been removed. Instead the user can compile the code with only the required modules, making for a smaller executable that is dynamically streamlined for their intended purpose as well as minimising compile time issues that may have slowed down users of the older code. Another advantage of a modular approach is that developers can add modules, and can avoid going through other subroutines in the suite to find where the new algorithm can sit comfortably and access the variables it requires during run time when developing code.
This modular approach allows for core modules to be treated as "black boxes" that can be put together like Lego bricks, simplifying development further. There is no need now for branching in run-time via “goto” or “if” statements, and adding new functionality to the suite simply requires writing a new module that can be "plugged in" to the code.
The release of more advanced versions of Fortran spurred an effort to rewrite the FIREBALL code from scratch, as part of the work for this thesis. The revised structure simplifies the development of future features in the suite and takes advantage of newer Fortran features. The project has been called “Lightning” after its speed increases and a more intuitive modular code base.
The initial and most obvious benefit is that the code has been greatly simplified. The previous (now called FIREBALL ’96) package continued to be added to and enhanced over the years. This resulted in a code that is up to date but had become quite monolithic, and difficult to develop further. The large code base required by, for example, the FIREBALL suite, can sometimes run into difficulty in compilation due to dependencies on libraries as well as processor and compiler issues. Such issues are difficult to debug in legacy code such as the ’96 package. Many improvements have been integrated into Lightning and the new Create (the program used to create the pre-computed tables discussed above) which achieves these goals admirably. Possibly the most effective aspect, is that the ’96 package used many "if" statements. "If" statements allow for the code to set up different calculations based on some state of the input. This allows for extensive user control over code operations by using input files (in plaintext) containing on and off tags. However, it makes the code difficult to follow as well as slowing down the overall runtime of the package. By employing modular programming, this is no longer an issue.
In the new package, output files are now "format free" making them easier to use outside of the package itself. There are now index files showing the user exactly which data is where on output.
Finally, extensive testing was carried out to optimise the methods employed in Create for carrying out integrations. The current implementation of these algorithms has proved itself to be both efficient in terms of computing resources and consistent in delivering accurate results.
The suite’s overall performance and operability have thus been enhanced for both users and developers. The complexity of the code has been reduced. This rewrite of FIREBALL has proved to be challenging, its results will show that it was worthwhile. Preliminary in-house testing indicates a significant increase in speed arising from the new optimisations within Fortran.
As a member of the team developing the Lightning package, this project entailed rewriting the entire Create suite as well as modules of FIREBALL-Lightning. This has included the evaluation of adaptive numerical integration techniques[47, 48] into FIREBALL and has allowed for the further development of more advanced features in FIREBALL.
This ongoing effort finally makes it possible to create a module for FIREBALL which will allow us to calculate an absorbance spectrum from the DFT results - a vital step towards the goal of full optical property calculations of solids. This will allow for the simulation of experimental techniques such as reflectance anisotropy spectroscopy (RAS), which makes FIREBALL a very powerful tool in the future for the characterisation of nanoscale systems and the comparison with experimental results on a direct level. This has already facilitated the development of an exact exchange module within FIREBALL, which is a major accomplishment towards the calculation of optical spectra.
The overall improvements in FIREBALL-Lightning make it far more portable than before as well as minimising its run-time. The simplifications introduced by derived type variables and modular programming make development for the suite far simpler than before. Within Lightning we now have a new tool that is incredibly versatile for energy calculations and simulations.
This chapter discussed the many aspects of ab-initio simulations that are particular to the FIREBALL methodology. The implementation of the methods used is discussed in detail and the differences between FIREBALL and FIREBALL-Lightning are explained. In the following two chapters the program which generates the Fdata for FIREBALL is discussed and the development of that program for the FIREBALL-Lightning suite is explained in detail.
As discussed in Chapters 2 and 3, one of the advantages of the FIREBALL package is its overall speed. In MD simulations it has been shown to approach an Order(n) scalability. Speed is achieved because all the interactions between atoms are pre-compiled into a dataset for the calculation known as the Fdata. Within the Fdata are a number of files corresponding to each interaction, calculated over a regularly-spaced numerical grid for two-center and three-center interactions. When the FIREBALL program itself is run, these tables are read into memory and interpolation is carried out to find interaction values for specific distances between orbitals as part of the MD or other such calculation. This chapter and Chapter 5 deal with the Create program which computes the Fdata.
The Create program is essentially the workhorse of the FIREBALL package. As discussed in Chapter 3, the matrix elements Hij and Sij (the Hamiltonian and the overlap) go to exactly zero beyond some cutoff distance (rci + rcj). This means that there is a finite and acceptable range over which the integrals of interactions need be calculated. Because these interaction tables depend only on the atom type, the chosen rc values and the type of exchange correlation functional chosen, they need to be calculated only once for a set of atomic species. As long as the criteria don’t change the Fdata can be reused for a number of simulations.
This precompiled Fdata approach contrasts with methods that calculate interaction integrals as they are needed, which can be computationally inefficient due to calculation repetition. It also contrasts with methods that evaluate integrals once per geometry and then store them to disk during the self-consistent procedure, leading to the disadvantage of extensive disk I/O (slowing the calculation) and computational cost. By evaluating integrals before any MD is carried out, FIREBALL’s Fdata approach attempts to reduce the disadvantages inherent in both these approaches. Within the FIREBALL program the pre-compiled interaction integrals generated by Create are loaded into RAM for fast access. When integrals that require evaluation are needed, an interpolation is carried out between data points taken from the tables in the Fdata.
The generation of Fdata within the Create program lends itself to parallelisation by breaking down the number of integrals required by the interaction type and the species. Each of these integrals may be evaluated on various nodes in a cluster. Parallelisation is important in the generation of Fdata because the number of data files generated is Order(N3) dependent on the number of electrons, N. By spreading this load out to a number of processors, calculation time can be greatly reduced. This scalability by parallelisation and over a number of species is shown in figure 1.
Within Create a number of interactions must be calculated and a number of optional interactions can be calculated. The following is a brief summary of these.
For two-center (2C) interactions, we evaluate everything in cylindrical coordinates. This simplifies many aspects of the calculation as we can then place the two centers on the Z axis. Beyond the limit of the some of the cutoffs of the pseudoatomic wavefunction, the integral is exactly zero, this is shown pictorially in figure (4.2).
The overlap integral is required for most interactions and is described in any introductory text on the subject of quantum mechanics. In FIREBALL we use the pseudoatomic wavefunctions to evaluate this:
where ψμ and ψν are two pseudo-atomic wavefunctions.
2. VNL (Voltage Non-Local).
The Non-Local Potential is due to the pseudopotential of the non-local shells of atoms. Non-local pseudopotential terms are expressed as the overlap between a pseudoatomic wavefunction ψμ NL and atomic wavefunction ψν:
3. VNA (Voltage Neutral-Atom).
The Neutral-Atom Potential is due to the pseudopotential of local atoms on a two-center interaction. As described in Section 3.2, there are three cases for VNA:
and the “atom” case:
where v is the potential operator and all other symbols are as described earlier.
4. Dipole interactions.
The integral of the general two-center matrix elements are of the form <ψ1 | r | ψ2>. These are the dipole terms. The z-dipole is required for the SCF calculations:
where z is the z-component position operator on the pseudo-wavefunctions, and thus describes the electron density in this direction.
5. Coulomb potential.
The Coulomb potential is simply the short-range coulomb potential using spherical densities. It is described by:
6. Kinetic potential.
The Kinetic Energy interaction matrix elements are calculated for two atomic species. The electronic KE is calculated in k-space to avoid numerically solving the derivatives, saving computational time. They are represented by:
where all variables have their standard meaning.
7. Average Density.
Charge density is required for both implementations of the exchange-correlation interaction within FIREBALL. This is explained in Section 3.2.2. Again we have three cases, those being ontop L/R and Atom.
where is the position operator.
8. Spherically-Averaged Density.
Spherically-averaged electron density is the same as that in 7 above, but with spherically-weighted wave functions rather than atomic wavefunctions. This is also needed for XC calculations.
9. Spherical Overlap.
Overlap using the spherical wavefunctions approximation.
Double-precision Extended Hubbard XC interactions, given by:
Where ρα is the electron density of one spin orientation, and ρβ is that of the other spin.
11. Exchange Correlation.
As overviewed in the previous chapter, the Horsfield XC interactions are given by:
where all symbols have their standard meanings.
Three-center (3C) interactions in Create are evaluated using a bond-charge scheme. The three centers form a plane (the π - σ plane). The σ axis connects two centers, where the orbitals μ and ν reside, with the origin assigned to the midpoint between them. The position of the third center is defined relative to this origin. This is depicted in figure 3.
As with 2C interactions, all interactions in 3C go to zero when a certain limit of distance is reached. In 3C, this is when dbc > rca + rcb. We again build Fdata files based on iteratively moving the atoms in figure 4.3 labeled a and b apart from dbc = 0 to dbc = rca + rcb and d = 0 to d = rcc. This is done for five angles of θ, which is sufficient for the interpolation algorithms in the FIREBALL program. The five angles for θ are chosen such that cos(θ) = 0, 1/√3, -1/√3, √(3/5) and -√(3/5). This results in five equations with five unknowns which are solved and the process is repeated over the range of values of d and dbc.
The three-center interactions that are needed for FIREBALL are:
1. Bond-Charge Neutral Atom (BCNA) which refers to the neutral and charged atom potentials. This is the 3C analog of the VNA in the 2C.
2. Den3C is the 3C density,as in the case of the 2C version. It is required for the OLSXC and SNXC XC options.
3. DenS3C is the 3C density in the spherical approximation, again, a 3C expansion of above.
4. XC3C is the Horsfield XC matrix terms in 3C
Prior to the improvements made as part of this work, the implementation of the Create program used input based on three major input files as well as an information file for each species in the study, giving four file types in total. The three major files are switch.input, theory.input and create.input. For the silicon carbide and indium on silicon studies carried out for this work, the old Create was used. For this reason, an overview is included here.
The first input file, switch.input contains a number of input toggles. It is a list of ‘1’ or ‘0’ for each of the interactions described in the section above. More precisely, it breaks down each of the inputs above into each of their constituent interactions (for example, vna ontop L =1, vna ontop R = 1 and vna atom =1), in total some 20 separate toggles.
The second, theory.input specifies, again by a number of ‘0’ or ‘1’ values, toggles according to which theories the user wishes to specify. DOGS can be toggled here, as well as HARRIS, SNXC, OLSXC. These are all described in the theory in Chapters 2 and 3. A sample of the theory.input is shown in figure 4.
The theory.input file is similar in format, but considerably shorter than switch.input. In this file, itest can in essence always be set to 0. Setting it to 1 means: “ignore the rest of switch.input” and is only used when debugging the code. The iharris, idogs, ihubbard and ispin tags turn on or off the Harris integrals, DOGS integrals, extended-hubbard integrals and spin-density XC interactions, respectively. The ixc_opt flag chooses between using SNXC or McWEDA (OLSXC) as described previously. The ioomethod toggle is currently redundant and the igauss3C option toggles on or off Gaussian fitting approximations to the three center interactions, reducing accuracy, but also greatly reducing run-time. This is usually done to check input parameters before generating a full Fdata set.
The third input file is the create.input. It lists the molecule files of each of the species for which we are calculating the Fdata.
The remaining input files are those referenced in create.input. The naming convention of these files is generally: ATOMICSYMBOL.PPTYPE.input where PPTYPE is the pseudopotential theory used, for example, the “3t” in the create.input above implies LDA (the ‘3’) within the Troullier-Martin PP (the ‘t’). For example, the Oxygen input file, (O.3t.input) is shown in figure 6.
The files referred to in this are the output from the pseudopotential generator and the Begin program, which is also part of FIREBALL. The Begin program generates the wavefunction, neutral atom, and non-neutral atom files, given by the extension .wf*, .na0 and .na*, respectively.
The Fdata files, the output of Create, have a naming convention based on the interaction type, species involved, shells of those species and in the case of 3C, an integer value (01 to 05) corresponding to the value of θ that the file corresponds to.
A flow chart of the old Create code is in Appendix V. The following chapter outlines the development of the new Create code; all interactions outlined above are implemented in the newer Create code, as well as optimisations of the algorithms involved in the new code.
This chapter has introduced the Create program and its methodology within the FIREBALL suite. The next chapter shows how this has been significantly enhanced during the course of this work for FIREBALL-Lightning. In Chapter 5 we show the enhancements achieved by rewriting the code base from the beginning and Chapters 6, 7 and 8 deal with the work carried out using this and the older FIREBALL code.
With any package like FIREBALL there comes a point when the code base needs to be rewritten, for further development to be possible. During the course of the research for this thesis, the FIREBALL suite was rewritten as FIREBALL-Lightning, which is a faster implementation of the code. This chapter discusses the rewriting of the Create program for FIREBALL-Lightning, the final result being the Create-Lightning program.
Much of the FIREBALL and Create code dated back to the FORTRAN77 era. Since the advent of Fortran90/95, many features have become obsolescent, meaning that they were marked for deletion in future Fortran versions.
Many of the advantages of newer versions of Fortran were outlined and discussed in the previous chapter. How these new features apply to the Create code specifically is next discussed here. For clarity, in this chapter the older Create program is referred to as Create-2006 and the new program is called Create-Lightning.
Firstly, there is a new automatic optimisation of loops which is very useful in Create-Lightning, due to the sheer number of iterations and integrations that must be carried out to generate an Fdata set. Where possible, the compiler can “unroll” loops, which makes them run faster by turning them into command lists rather than an actual loop. Hard-coded loops that have been unrolled are more efficient as there are less branches for the computer to execute in run-time.
The previously-described array operations also remove the need to specifically code in loops to work through an array. As mentioned before, this also means that at compile-time, architecture-specific optimisations4 can be implemented. This also applies to structs, which make data handling in Create-Lightning far more intuitive.
Secondly, for each interaction described in Sections 4.2.1 and 4.2.2, we now have a specific plugin module. This eliminates the toggle files that are discussed in Section 4.3.1 and required for Create-2006, simplifying both the usability and extensibility of the code.
Thirdly, a feature which was tested, but eventually not used, from newer Fortran versions is recursion. Recursion in this context is the ability of a function or subroutine to call itself. It is a feature first added in Fortran90. It allows for adaptive quadrature to be implemented because it allows for a function to call another incidence of itself until some parameter is met.
The last addition to the Fortran90 specification that should be mentioned at this point is the concept of interfaces. Previously, in FORTRAN77, it was only possible to pass number variables and character strings between functions and subroutines. By defining an interface, a function or subroutine can be passed as an argument to a different function or a subroutine. As a tool, this is one of the most powerful features of the new Create-Lightning program. It means that instead of having reused code in every module for an interaction, the core code can contain subroutines to carry out, for example, the grid iterations and integrations needed for all interactions; then, only the actual mathematical function being evaluated need be passed in as a piece of code, and all other tasks are completed around that function.
These newer Fortran features are invaluable in making code run more efficiently. Where possible, they are used to their fullest potential in the rewrite of FIREBALL. In essence, the new code is far easier to read than the monolithic older implementation. It is clearer and simpler to add new features as modules, and the code runs faster than before, even on a single processor machine.
The rewrite of FIREBALL is outlined in the previous chapter. This section deals more specifically with the Create program. Where possible we made use of the newer features in Fortran95 and later releases, such as array functions, derived types, and modularisation.
The entire Create code was rewritten during this work. As well as the modularisation of interactions mentioned in Section 5.1, other changes in the new Create-Lightning code are that the Fdata files now have no headers -instead a set of index files is generated. There are a number of reasons for this. First, it makes it easier to make Create-Lightning append an already existent Fdata set, should the user wish to add more interactions or species. Second, by using index files, the user can see far more easily what has been compiled in the Fdata, thus not relying entirely on a naming convention adopted for the Fdata directories as was previously done. In run tests the new Create-Lightning code completed a generation of a full Fdata in about 85% of the time taken by older versions when run on a single processor.
For each of the modules written to calculate the interactions described in Sections 4.2.1 and 4.2.2, the output was compared directly to the output of Create-2006 to ensure that they both corresponded. Furthermore, an Fdata from Create-2006 was modified by removing the headers from files to allow it to be read by FIREBALL-Lightning, and results from new Create-Lightning/ FIREBALL-Lightning were compared with Create-2006/FIREBALL-Lightning and Create-2006/FIREBALL-2006. This procedure was required to ensure that all new modules were working as expected. At time of writing, no uncorrected errors have been found. Each module contains only what is unique to the interaction, and calls core subroutines to carry out the rest. This minimises the code base, which is an advantage in debugging and simplifying the usability of the code.
A Fortran program was written to run through all values within the output files directly from both Create-Lightning and Create-2006. The output files within each Fdata were compared. This was possible because all Fdata files are on a real-space grid and should be equivalent. If the calculated result at any point for any interaction differed by some small amount, such as 0.001%, a flag was raised. This assisted in finding bugs and confirming that the new Create-Lightning program gave the same results, shown to be correct in a number of previous publications, as the Create-2006 program.
As part of the effort to speed up the code, an attempt was made to add a new integration routine, based on adaptive Simpson’s quadrature.
Simpson’s rule is standard in any basic mathematics text book. For the most part, within Create-2006, we had employed a Simpson’s rule numerical integrator with the interval set at d/107, where d = rc1 + rc2. The dividend of 107 was chosen by the authors of the original Create-2006 program and has survived by legacy, having been found by trial-and-error to be a balance between computational accuracy and computational speed. Throughout the history of the Create-2006 program, this appeared sufficient. However, with the Create-Lightning, it was felt that the time had come to attempt to update this to something more modern and efficient, such as an adaptive Simpson’s implementation.
The adaptive Simpson’s implementation is only possible within Fortran90/95 and more recent Fortran versions because of their ability to pass a function to a function. Thus a routine can be sent to a function called “adap_simpson”, for example, along with the limits of integration and tolerance arguments.
The adaptive Simpson’s method is based on estimating the error arising from calculating a definite integral using Simpson’s rule. If the error exceeds a user-specified tolerance, the algorithm calls for halving the interval of integration and applying the adaptive Simpson’s method to each subinterval in a recursive manner. This continues until:
where a and b are the ends of an interval with midpoint c, S() is the Simpson’s rule estimate of that interval and ϵ is the lowest non-zero number that the computer architecture can resolve, which is dependent on the compile parameters. The 15 ensures that estimates obtained are exact for polynomials of degree 5 or less. (This is shown in ref. []). The value of ϵ in single precision is 1.1920929E-07, in double precision is 2.220446049250313E-016 and in quadruple precision is 1.925929944387235853055977942584927E-0034. This corresponds to the real numbers being 4, 8 or 16 bit, respectively, on the Intel X86 architecture used in the Lewis Group’s computational cluster.
The “adaptive” part of the algorithm works by looking at the result of 3 evaluations of an area. The first is to numerically evaluate the integrand between two points. The second and third evaluate half that interval each. If the result of the first is greater than a specific tolerance of the sum of the second and third, the program calls itself with tighter limits of integration, splitting the intervals into even smaller sections. Once tolerance is reached these smaller intervals can be summed to find the result.
There are two main advantages to this method over a standard Simpson’s calculation. Firstly, because it is “hunting” until it reaches a specific tolerance, for well-behaved functions there is less computation required to evaluate a function. Secondly, it is more exact than a standard Simpson’s rule implementation as it is not limited by the set interval parameters inherent in standard Simpson’s rule.
As a first step, an adaptive Simpson’s algorithm was implemented within a stand-alone program and tested extensively. For simple trigonometric functions the adaptive Simpson’s routines worked smoothly and effectively. The subroutine produced results that were within 0.01% of analytically derived results up to polynomials of order 6 without any issue. After that, results were within 1%, which is acceptable as adaptive Simpson’s algorithm is no longer accurate at this order of polynomial. Once we moved to polynomial expressions of order > 8, there was appreciable discrepancy between exact solutions and the function’s solution. (A table of these results can be found in Appendix I.)
The algorithm was then implemented within the new Create-Lightning program for all integrations. It was also applied to the Exact Exchange that was implemented in Create-Lightning, which will be discussed later. Within the new modular framework, evaluating different integral techniques is not a difficult matter as we can plug in a new algorithm easily.
On testing the Simpson’s method, there appeared to be an underlying issue in results which is still undiagnosed. This issue became obvious from the appearance of “numerical blips” in output data. It appears that calculations run cleanly, and there are no obvious errors in the code. Checks were done to ensure that there was not a race condition or a divide by zero error. These attempts were not successful in diagnosing the problem.
The graphs in figures 5.3 - 5.6 are from one interaction within Create, the density ontop L case, in which the issues were first detected when testing the new Create-Lightning code on this module. It is not isolated to this case. However, this case was the first appearance of the problem and it was quite well pronounced. Figure 5.3 shows this issue clearly.
As can be seen in the Adaptive Simpson’s Quadrature case, there is a noticeable “blip” at both points 4 and 20. This corresponds to two-center atomic distances of 0.7935Å(3*0.2645Å) and 5.0255Å (19*0.2645Å), respectively. This is again obvious in figure 5.4.
As part of the debugging to track down the issue that caused the discrepancy between the adaptive Simpson’s and the standard Simpson’s algorithms, a verbose version of the module and Create-2006 was added for the interaction that figures 9 and 10 produced. The debugging version showed that for some of the points on this graph, the Adaptive Simpson’s Quadrature only required 3 or less iterations. It is, however, a recursive algorithm, so it sometimes required over 300 iterations. By comparison, standard Simpson’s requires a hard-coded 213 iterations ((107 * 2) -1).
The outputs shown in figures 9 and 10 were carried out using single precision. The functions can be compiled with the “-r8” and “-r16” flag for double and quadruple precision. By using higher precision, the discrepancy between any two values can be smaller and it was postulated that this may be the cause of the errors in figures 9 and 10, however the results from the adaptive quadrature were very similar to what is seen in these graphs.
In another attempt to diagnose the errors seen in figures 9 - 10 it was assumed the adaptive routine was correct and that the standard routine “smoothed over” some inconsistency in the test input data set, given that the standard Simpson’s algorithm is fairly coarse. Extensive graphing of the interpolated wavefunctions, (i.e., the input to the adaptive Simpson’s subroutine) yielded no such anomaly. Neither did a verbose output of all input parameters during runtime.
A number of attempts to find the source of the errors are outlined as follows:
Experiment 1. Use double precision.
Experiment 2. Create a “Simpsons” subroutine, implemented in the same manner as Adaptive Simpson’s. This was to ensure that the error was not a message-passing problem of some sort.
Experiment 3. Increase the number of hard-coded iterations that standard-Simpson’s uses. This was an attempt to recreate the error produced by adaptive Simpson’s.
Experiment 4. Use quadruple precision. Experiment 5. Invoke a different quadrature subroutine- this was taken from the scientific library QuadPack.
Results of these tests are as follows:
Result 1. Both sets of code, the non-Simpson’s and the Adaptive-Simpson’s were compiled with the “-r8” argument using the ifort compiler. This flag tells the compiler to make all variables double precision.
A comparison of both non-adaptive and adaptive algorithms with this flag in compiling the code showed no change in the non-adaptive output, which is not surprising. In the adaptive case, there was a slight smoothing of the output, but the “blip” shown in figure 5.3 at point 20 remained.
Result 2. In an attempt to see if there was an issue with variables passing to the adaptive Simpson’s module, an additional subroutine was added to that module. This was modelled on the “adap_simpson” subroutine in terms of input, but only carried out standard Simpson’s integration.
Results from this were identical to the non-adaptive Simpson’s results seen previously. This rules out any message-passing errors. This is shown in figure 11.
Result 3. Further work involved attempting to increase the number of grid points used by the non-adaptive Simpson’s method. The grid size was increased from 107 points to 187 points, then to 424. These values were chosen randomly.
Despite marked increase in time for the computation to run, the results were, again, identical to what was previously seen in the standard non-adaptive Simpson’s runs.
Result 4. The three studies outlined above, however, do indicate that there is an issue with the adaptive Simpson’s algorithm as implemented. Working from the results of Experiment 1, the code was compiled with the “-r16” flag, thus making all integers 16 bit numbers. Again, standard Simpson’s showed no change in results, and again adaptive Simpson’s showed a blip in the same place as in Experiment 1.
This still indicated an error within the adaptive Simpson’s subroutine, and careful analysis of its working variables data was carried out. This is outlined in Appendix I, but essentially, some of the calculation variables can reach a very high or very low value at times, and in double- or quadruple- precision these values sometimes differ greatly from that in the single-precision mode. This indicates that the limits of the computer are being reached. However the cause of this remains elusive. Due to the tolerance setting which prevents having to deal with any number below 15ε, the criteria shouldn’t be met.
Result 5. The use of a standard-library quadrature implementation in Create-Lightning to carry out the numerical integrations was considered. This was carried out by employing the QUADPACK library (www.netlib.org/quadpack/) and then using it in place of the “simpson” or “adap-simpson” subroutines. Quadpack is part of the Netlib Library which also contains SCALAPACK, which is currently used in FIREBALL. This was now easy to implement in Create-Lightning because of the modular setup.
The applicable subroutines are QAG or QAGS in this case. QAG is a simple globally adaptive integrator using the strategy of Aind, ref []. It is possible to choose between 6 pairs of Gauss-Kronrod(GK) quadrature formulae for the rule evaluation component. QAGP is a polynomial-optimised GK integrator, and QNP is a Non-Adaptive quadrature subroutine.
There are two user control variables for all of these subroutines, relative error and requested absolute error. The subroutine outputs the estimated absolute error during runtime.
Originally issues arose from the implementation of QAG. The issue is that at certain levels of relative and absolute error, the output flags show a number of errors. These were mainly error codes 2 and 3, which correspond to the output statements “the occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved” and “extremely bad integrand behavior occurs at some points of the integration interval”, respectively.
Results are shown in figs (5.6a-f). Any adaptive subroutine that was applied with an acceptable error resulted in integrations that were not useable. Both QAG and QAGS were tested, which are general integrator routines, QAGS being optimised for sinusoidal functions. This data strongly suggests that the accuracy of the desired output is dependent solely on absolute error, and an absolute error of only 0.1 is not within acceptable margins.
This trend appears to continue, as absolute error approaches or goes beyond something useful- for example 10-3. The output then becomes unusable. A number of reasons have been proposed for the emergence of these errors. One is that “catastrophic cancellation” errors may be occurring. This is when two close values are subtracted from one another and the answer is incorrect for numerical reasons. This is hard to accept given that such well established libraries as the QUADPACK routines were being used.
With the issue unresolved, the standard Simpson’s integration algorithm continues to be used in the Create-Lightning program rather than the quadrature approach. Implementation of a better quadrature itself is a very trivial matter due to the modular nature of the new Create-Lightning. Despite the newer algorithm not being used, it is felt that due to the very nature of the lower value maxima of the wavefunctions, the limiting factor is the computer hardware and not the choice of integration algorithm.
An underlying aim of this work was the potential development of a system for the calculation of optical spectra. As part of this proposed development, exact exchange (EXX) would be needed and was therefore implemented in Create-Lightning. Use of the generated Fdata is currently being implemented in the FIREBALL program.
To properly evaluate optical properties, excited states need to be describable. It is well known that the LDA does not describe energy gaps well. One method for evaluating these gaps in a truer form is to use the so-called exact-exchange density functional theory which allows the determination of the exact local Kohn-Sham potential. It has been studied in much detail for semiconductors, for example in refs [[52, 53, 54]].
Implementation in Create-Lightning was done by making use of the new modular system. This had the added advantage of allowing for testing of the ease of implementation of new modules as well as making the new module required.
For Vex being the exchange potential, we can write the exact exchange potential as:
where e is the electronic charge, r and r’ are the center positions of wavefunctions Ψi and Ψj, respectively.
As before, we are interested in the expectation value between two atomic orbitals, which results in a four-center expression, within the LCAO formulation of (see Appendix II):
where Vexν indicates the potential operating on the atomic orbital. All other symbols are as before. This is a four center integral, and we simplify as we did before in Chapter 4 for other interactions with one center, two center and three-center cases. The full four-center integral is not required, as its relative contribution is small.
The final solution to the exact exchange algorithms is shown in equations 5.4 - 5.7, with a derivation in Appendix II.
For the one-center case we get:
which is arrived at from equation 5.3 by simply letting r1 = r2 = r3 = r4 and α = β
For the two center cases, as in the general theory chapter, we have an atom, ontop L and ontop R case, given by:
Atom case, r1 = r2
Ontop L case, r3 = r1
Ontop R case, r3 = r2
Unlike previous exchange approximations explained in Section 2.2.5, EXX is exact. We deal with the Coulomb kernel by replacing the 1/|r - r’| term with a double sum of associated Legendre polynomials multiplied by cosine terms, using the addition theorem for spherical harmonics. We simplify these in code by analytically using Clebsch Gordon and Legendre Polynomials on the angular part of the wavefunction.
The exchange potential isn’t a directly physical phenomenon and it is difficult to test the output of this module directly. It can be seen by inspection that the output behaves well. There are smooth curves due to the Coulomb kernel. Testing of the output of this module has proven difficult. Currently the only way to test the output is to compare it with code that carries out the same interactions. This was done by using the implementation of Daniel Jensen, from Brigham Young University. Comparison graphs of these outputs can be seen in figure 5.7. This sample is from the two-center case where there are two orbitals on each center.
EXX is an invaluable addition to the FIREBALL suite for calculating optical spectra. The band structure generated by exact exchange is a far better approximation of the actual band structure. With LDA and GGA, it has been shown that the band energies are too close together with respect to an experimentally designed result, whereas within the EXX, band energies are correct. It is for this reason that the generation of optical spectra, due to electron excitation, is required.15 The EXX module has been completed for Create-Lightning as a plugin. Its runtime is significantly longer than the other interactions as it scales by O(N4) with the number of electrons N in the system.
The EXX can be combined in FIREBALL with the LDA correlation terms to yield a quantitatively more accurate result than LDA exchange-correlation.
Within the FIREBALL program itself, it is crucial that linked lists of atoms and their neighbours are maintained. For each time-step, the forces incident on each atom are a function of the surrounding atoms that fall within a distance d (= rcatom + rcneighbour). Within the new framework of FIREBALL-Lightning, this required a rewrite of the “neighbors” module.
This was the only module written by the author for FIREBALL-Lightning, and was carried out at an early stage in this work as an exercise. It generates tables of neighbours of each atom, and a reciprocal table of this data. The reciprocal table is required for bookkeeping so as to avoid double-sum interactions.
A rewrite of the Create code was a considerable undertaking while progressing the work for this thesis. Section 5.3 describes how the code takes advantage of modularisation and other features of Fortran95, and the advantages this brings for further code development. The speed gains of the code are shown in figures 5.8 and 5.9, which show the time of calculation for systems of increasing complexity and for an increasing number of nodes in parallel processing.
Up to this point in the thesis, the background and implementation of the FIREBALL suite have been discussed. The theory that underpins the code, and then the code itself have been described. The next three chapters deal with this code as applied to real systems which can be compared with experimental data from the literature. In Chapter 6, the use of MD to explore a semiconductor-metal transition due to deposition is explored. Chapter 7 describes a temperature-induced phase transition and its analysis with FIREBALL. Chapter 8 discusses the use of FIREBALL in evaluating the optical properties of a system.
This chapter and Chapter 7 discuss work carried out on silicon carbide surfaces. Section 1 of this chapter serves as an introduction to silicon carbide in general for both chapters. Section 2 functions specifically as the introduction to the potassium surfaces on silicon carbide. Sections 3 and 4 discuss the experimental work carried out to understand this surface better.
Silicon carbide (SiC) as a general material has a number of uses. It is used as a gemstone, known as Moissanite, which is named after the Nobel Laureate who first discovered it in nature in 1895. It has mechanical uses in applications such as abrasives, saw-blade coatings, bulletproof vests, as the reflective coating on special application mirrors (such as in the Herschel Space Observatory) and it was also used to make the first blue LEDs.
In electronic applications, SiC exhibits a number of interesting properties and has been the subject of a large number of studies. For example, due to its relative ease of manufacture, it may be possible to tailor-make devices with specific properties, e.g. controlling the band gap size. It has been shown that it is very stable, and can resist degradation from radiation or other harsh environments. It has been shown that the switching rate of an SiC transistor may be capable of far exceeding that of Silicon technologies. It is a wide bandgap semiconductor and is, unusually, a IV-IV non-elemental semiconductor. As a result it is interesting since it can be studied as something in between the elemental semiconductors and the III-V’s. The properties of SiC as a semiconductor have been explored in depth, due to the uniqueness of its chemistry which shows some exceptional behaviour.[60, 61, 62, 63, 64]
SiC can exist in a number of energetically degenerate polytypes. Depending on the polytype, the bandgap has been measured from 2.38 to 3.26 eV. It also exhibits some very interesting behaviour, such as surface metallisation upon hydrogenisation of the 3C-SiC(100) surface and a Mott-Hubbard transition on the 3C-SiC(111) surface and equivalent hexagonal surfaces.
This chapter serves as an introduction to silicon carbide as a material and discusses the work carried out for this project on the potassium-induced reconstructions of the 3C-SiC (100) surface. The FIREBALL package is an ideal tool for studying problems of this nature, since, owing to its computational efficiency, many hundreds of possible surfaces can be studied. Rather than PW codes which require far more computational resources for similar simulations, FIREBALL can reuse its Fdata for each surface topology studied, greatly reducing the computational cost of such work. In Chapter 7 the research carried out on the clean 3C-SiC(111) and 6H-SiC(0001) surfaces is discussed, which exhibits a “soft phonon” interaction. FIREBALL is also well suited for soft phonon studies because very long runs, i.e. a very large number of time steps in MD, can be carried out with great efficiency.
Silicon Carbide can exist in a number of different physical forms. There is an amorphous phase, as well as at least 200 differing crystal polytypes.[68, 10] The classification of these polytypes by their general topology is simple. By defining three different forms of bilayer, A, B and C, these polytypes can be uniquely identified.
In this classification system, types are described by the number of bilayers in their primitive repeating unit. A letter then describes the symmetry of the unit, for example, 3C-SiC has three layers in its repeating unit cell, and exhibits cubic symmetry, whereas 6H-SiC has six layers in its repeating cell and is hexagonal. The other symmetry notation is “R” for rhombohedral.
The descriptions of these bilayers a, b, and c are quite intuitive. In their simplest form, each of their units consist of three atoms and the bonds between them. For example bilayer a is depicted by the image in figure (6.1a). Bilayer b, in figure(6.1b), is nearly identical to a, except for the fact that it is a translation of bilayer a. This subtle difference is more apparent if both a and b are shown together, as in figure (6.1c). A polytype of SiC, known as 2H-SiC is described by these two building blocks alone. It is also known as pure hexagonal SiC. This is the wurtzite structure.
The final building block that is needed for the description of polytypes is bilayer c. The c bilayer is generally similar to a and b. However it is rotated in the plane by 60 degrees with respect to a or b, as can be seen in figure (6.2a). Figure (6.2b) shows the three blocks together, to try to clarify the cataloging method used. This notation system is known as Ramsdell notation. The 3C structure, seen in Fig(6.3a) is the well-known zincblende structure.
Of interest in these studies are the 3C-SiC and the 6H-SiC polytypes, as shown in figure 6.3. There are a number of reasons to focus on these two polytypes specifically. Experimentally the 3C and the 6H polytypes are quite accessible, which has fueled a lot of work in the literature. Due to the interesting electronic properties and the accessibility to experimentalists, this work examines potassium deposition on the (100) surface and the Mott-Hubbard transition reported on the (111) surface.
In 2001, Derycke et al. published a surprising result from experiments on the (100) surface of 3C-SiC. Using STM and ultraviolet photoemission spectroscopy they discovered that by adsorbing hydrogen, the semiconducting surface transitions to metallic. This was surprising because previously it had been shown that hydrogenation of any previous semiconductor surface studied would passivate all dangling bonds, rendering the system an insulator. Such behaviour currently appears to be unique to the 3C-SiC(100) surface.
These results are surprising and yet have no definitive theoretical explanation. The problem has been picked up by a number of theorists. Using ab-initio methods Derycke et al.’s results in the electrical properties have been confirmed[69, 70, 71, 72] but there is still a debate over the proposed topology of the system. Theoretical results suggest that the surface hydrogen resides between atoms in the third layer of the surface in some sort of 2-electron atom bond. This is in contrast to experimental work using STM, STS and synchrotron photoemission spectroscopy suggesting that the hydrogen is directly on the surface. The theoretical studies therefore indicate that the metal transition happens on saturation of the Si-rich 3C-SiC surface, whereas unforced hydrogenation does not seem to have the same effect.
These results prompted Derycke to look further into the phenomenon of metallisation of the surface. He postulated that the addition of a group-I metal might also show some interesting results. In contrast to hydrogen, group-I metals have been shown previously  to make a semiconducting surface metallic. Alkali metals on semiconducting surfaces have been studied extensively. These studies show strong evidence for the idea that the weakly-bound s electron transfers to the substrate, resulting in metallic behaviour at saturation. As a result, it is interesting to explore these metal-semiconductor interfaces.
In this work, the surface reconstructions suggested by experimentalists were studied computationally. The results were compared with an in-depth study of other possible surface reconstructions and an evaluation of the electronic properties of these variations. Uniquely to the 3C-SiC(100) surface, it had previously been shown that the addition of potassium does not modify the silicon-rich c(4 × 2) reconstruction seen on the clean surface. This reconstruction has been dubbed the alternate up-down dimer model (AUDD). Other adsorbates have previously been shown to destroy the 4 × 2 AUDD array.
The structure of the clean c(4 × 2) surface was under contention for some time, however it is now accepted that the AUDD model is the correct one. It consists of alternate up and down dimers in a regular pattern. A number of alternative structures for this system have been proposed, for example the missing row asymmetric dimer (MRAD) model as well as others. Recent results, both computational and experimental agree on the AUDD model as being the most favourable overall.
Figure 6.4 highlights the position of "pedestal" sites on this surface. These are positions identified experimentally where additive potassium atoms appear to reside, based on results of Derycke et al. using STM, STS, LEED and synchrotron radiation-based core level photoemission spectroscopy (CLPS).
Experimentally, the added potassium preferentially fills these pedestal sites. Experimental results suggest that these sites are initially grouped together in pairs, forming a 2 × 3 pattern with two thirds of the sites occupied and one third unoccupied in a regular pattern, as shown in figure 6.4(b).[61, 62] As mentioned earlier the addition of potassium does not appear to affect the underlying surface and the overall surface remains semiconducting at low coverages.6 Other adsorbates have previously been shown to considerably alter the underlying structure. Fig(6.4b) shows the general schematic of this reconstruction.
With additional deposition, this becomes a 2 × 1 surface with all pedestal sites filled (see figure 6.4c) that is metallic. The experimental indication is that the potassium atoms "short circuit" the semiconducting nature by creating a link between all surface atoms at saturation. The underlying surface appears to remain unchanged, which is almost unique for a metallic overlayer on a semiconductor.
Derycke et al. report that from LEED results, with deposition of potassium, the surface undergoes a change from c(4 × 2) to (2 × 3) to (2 × 1) at saturation. This is similar to the behaviour observed with the deposition of hydrogen - the c(4 × 2) SiC surface transforms into a (2 × 1) reconstruction. SXRD results show that the electronic properties of the 3C- SiC surface are considerably modified by the addition of potassium. At lower coverages as exhibited by the (2 × 3) reconstruction, the surface remains semiconducting, whereas with the saturated (2 × 1) reconstruction, it is metallic.
There is currently no experimentally-derived STM image of the potassium induced surface available. However, published STS results which can be considered to represent a local density of states at the surface, show the bandgap closing at θ = 1ML and acting as a metal.
That the underlying structure does not become modified by the addition of potassium atoms is unique. Ref [] reports a charge transfer from the potassium surface to the SiC underlayer, indicative of physisorption. The reported potassium-potassium distance on the surface is 3.08Å. Given that experimental results for potassium-potassium distances in a potassium molecule (K2) are between 3.84Å and 4.22Å with an average of 4.04Å, this would indicate an overlapping of the s and p orbitals of the potassium due to their proximity. It is possible that the potassium is lying on the surface and forced into a configuration that allows for electron transfer and conduction as a result of this proximity.
Experimentally, a smaller charge transfer has been indicated between the adsorbed potassium atoms and the surface than that between hydrogen and the same surface. This is experimentally evaluated by measuring the work function of the potassium, which decreases to about -3.2eV (from 2.29eV for elemental), which is very similar to that of potassium on the Si(100) surface.
These phenomena were studied within FIREBALL, and it was found that the surface potassium atoms do not form dimers as previously proposed. Rather, a new structure for the low-energy reconstruction in the θ = 1 monolayer coverage has been found, and it is shown that this structure does not alter significantly with large temperature changes. Lower potassium coverage on the surface tends towards the formation of zigzag chains and the potassium reconstruction is not dissimilar to bulk-like conditions.
Using FIREBALL and an optimised basis set, a number of calculations to look at the potassium induced surfaces in depth were set up. A minimised basis set of Fdata using hydrogen, carbon, silicon and potassium in the LDA was generated and used.
For hydrogen, a single s shell, with 1 electron, and a wavefunction cutoff of 4.50Å. was used. An excited shell of zero occupation, cutoff of 4.50Å and secondary quantum number of 0 was defined also. This allows for charge transfer to the hydrogen when simulating bulk.
The carbon was defined by two valence shells- the s and the p - with cutoffs of 4.20 and 4.40Å, respectively. Using the method defined in ref [], an occupation of 1 and 3 on these shells, respectively, was chosen. The silicon consisted of an s and p shell of cutoff 6.00Å each with an occupancy of 1 and 3. Finally the potassium consisted of a p and an s shell, of occupancy 6 and 1, with cutoffs of 6.00Å.
These parameters represent an optimised sp3 basis for the Si, C and K and a double ss* basis for the hydrogen. The method of optimisation is discussed in more depth in ref.[]. The SiC supercell consisted of 240 atoms, which constitutes 4 layers of SiC with a (100) surface. The lattice vectors were set such that in-plane reflections were accounted for, setting up a simulation of an infinite plane. The vector in the third direction, out of plane, was set to 999Å, to avoid any interaction between reflections in this direction. The bottommost layer was terminated with hydrogen and the atoms were fixed in position to simulate bulk. This constitutes a surface size of twice the c(4 × 2) or eight times the surface size of (2 × 1). With reconstructions that required further analysis, this supercell was doubled for further studies so as to double the size of the surface.
Unless otherwise stated, for all experiments the system was allowed to run MD at 300K in 2000 timesteps of 0.25fs, then quenched over ~600 timesteps. The resultant structure was then analysed. This process is known as simulated annealing. Many of these results were later confirmed by MD over longer time periods and quenching operations which yielded equivalent results.
A series of tests to confirm experimental parameters such as bulk modulus, were carried out and these results matched those of Trabada et al., for a 3C-SiC(100) surface with no adsorbed potassium. Tests were run to choose an optimal position upon which to place the potassium in the (100) direction, i.e. how far above the surface it could be placed. This also allowed us to see if the potassium has a preferred position under or in the plane of the surface. This meant that calculations could find equilibrium more quickly by the potassium being initially placed close to an optimum height above the surface. Figure 6.6 shows the supercell structure that was used.
Except for unphysical cases where the potassium layer contained more atoms than the surface could support and, for example, expelled the excess atoms, the potassium layer typically came to rest at ~2.2Å above the surface. This set up an ideal starting point for all subsequent calculations. A series of ~140 separate calculations was carried out to determine this starting height above the surface, varying displacements in the (100) direction above the surface and varying configurations and percentage coverages of potassium.
A number of calculations were set up to characterise the behaviour of two potassium atoms on the surface, representing a very low coverage. Experimentally, it is claimed that potassium initially forms dimers on the surface at coverages of < 2/3ML. Simulated annealing of these surfaces was carried out and a picture built up of how far each of the two potassium atoms come to rest from each other. The minimum potassium-potassium distances and total energy are graphed in figures 6.7 and 6.8.
In another set of studies the surface was loaded with potassium atoms one at a time, based on a set of loading parameters as outlined below. A potassium atom was added to a surface, and an annealing procedure was carried out, which was then repeated. This was done up to a coverage of 1ML, then in steps of 2 atoms up to 2ML. The potassium atoms were loaded based on one of four criteria, leading to four separate studies:
1) Placing the added potassium atoms as far away as possible from one another and, where possible, with a further row of pedestal sites between them.
2) Building up potassium atoms in the direction of trench growth, with dimer placement based on experimental results, while tending towards chains in this direction as the number of surface potassium atoms increases. This series also included slight modifications in this formulation - some included chains that were not translationally symmetric, as suggested in ref [].
3) As in 2) above, but with chains tending in the direction perpendicular to the chains in 2).
4) Random placement of dimers on the surface, orientated in various directions. A table and scatter plot of experimental parameters and total binding energy for all of these systems is in Appendix III. A chemical potential plot for the potassium was generated for the lowest energy structures found. The chemical potential of potassium is taken from calculations of a single potassium atom with vectors arranged so as it is in a BCC lattice. This value was calculated to be -435.83079eV. These are in the results section in figures (6.9) and (6.10).
In order to calculate the binding energy, the single-point energy of the supercell without any potassium, the energy of the potassium from above, and the total energy of the system in question, respectively, are used. The binding energy is calculated as:
where EB is the calculated binding energy, ESiCK is the energy of calculated SiC-K supercell, ESiC is the energy of calculated SiC with no added potassium, n is the number of potassium atoms added, and μK is the per-atom energy of potassium in the bulk, as in the previous paragraph.
where EB is now a function of μK, which is then plotted for the most favourable systems.
Alongside these studies, an examination of the experimentally-derived surface reconstructions for θ = 2/3ML and θ = 1ML was carried out. An adaptation of the θ = 2/3ML surface was also modeled, by shifting the dimers by one row in the direction parallel to the dimers.
As a final step the most likely reconstructions that were found for this surface, based on total energy calculations, were selected and the results rigorously retested against the experimental models. A larger supercell was used, to dispel any possible effects of the repeating supercell within the code, thus making the unit cell a square. MD simulations at 300K, 600K and 900K were also carried out. Higher temperature simulations were carried out for a number of reasons. Firstly, to see how stable these reconstructions were. Secondly behaviour of the surface at higher temperatures is indicative of the binding involved. And thirdly, higher temperature simulations allow us to see if any of the reconstructions may change completely given sufficient energy, as is the case with a temperature induced reversible phase transition.
For all annealed surfaces a density of states (DOS) and a surface band structure plot were generated. A metallic DOS for coverages of 1ML and a semiconducting DOS for any coverages below 1ML were expected.
The graphs in figures 6.7 and 6.8 show the results of the surface-dimer calculations outlined in 6.2.2. The distance between the two potassium
atoms on the surfaces is graphed agains the cohesive energy of the system, as described in equation 6.1. Figure 6.7 shows the chosen surfaces when simple quenching was carried out, i.e. the temperature is reduced in the simulation until a minimum is found. In figure 6.8 MD is carried out at 300K for 2000 0.25fs time steps before quenching, allowing the surface to reconstruct, this process is simulated annealing.
In figure 6.7, where the surfaces are simply quenched, the lowest energy topology that was found was the one labelled “j”. If ignored as an outlier, the graph in figure 6.7 follows a very pronounced curve, indicating that the potassium atoms do not interact well on this surface. The topology labelled “j” was a test based on the experimental result that dimers preferentially form on the surface in the direction that chains form in. In this case, “j” is in the perpendicular direction to the prescribed direction of chain formation, whereas “b” is in the direction shown experimentally. The result in “j” implies that at low coverages dimers do not form in this direction.
The similar energy between “b” and “c” is expected, since they are equivalent when the supercell approximation is considered. However, the direction of growth described by experimentalists forms dimers in ways more similar to “d”, which results in a structure that is unexpectedly of higher energy. This re-enforces the result with respect to surface “j” in the previous paragraph, that dimers, if they form at all, appear to prefer the perpendicular direction to those implied in the literature.
Figure 6.8 shows the same initial surfaces. However they have been subject to simulated annealing rather than simply being quenched. This process allows the atoms in the supercell to move for some time before cooling. Comparison of the images of the surface shows considerable reconstruction in the potassium between the quenched and the annealed surfaces. Note that the energy scale in figure 6.8 is lower than that in 6.7, which shows that all annealed structures are far more energetically stable than their quenched counterparts. As can be seen in figure 6.8 also, the relationship between potassium-potassium distance and the cohesive energy is still apparent.
A number of topologies based on differing conditions of placement were next studied, ranging from a potassium coverage of 0.083ML to 1.25ML, this corresponds to 1 potassium atom on the supercell surface through 15 potassium atoms on the surface. The scatter plot in figure 6.9 is a graph of the coverage Vs the cohesive energy calculated for each, after simulated annealing. In figure 6.9 the blue green X’s show the calculated cohesive energy for the LEED-derived surface structure from ref . As predicted experimentally, there is a minimum energy about the θ = 1ML coverage. There is a very large energy difference between ref []’s 1ML structure and the surface structures that were also tested, which is a compelling result and is explored in more detail in this section.
The scatter plot in figure 6.9 maps results from studies of coverages between 0ML and 1.25ML. At coverages of 1ML and above, extra potassium atoms were expelled from the surface. Exploring the large energy discrepancies between the experimental 1ML reconstruction and the results presented here is best done by limiting the number of reconstructions focussed on. The lowest energy reconstructions from figure 6.9 were selected and a chemical potential plot for potassium was generated, which can be seen in figure 6.10. There are a number of features of note in these two graphs. Firstly, the experimentally derived θ = 1ML coverage reconstruction is, from these results, far less favourable than some of the other reconstructions this study has identified for this surface. Secondly, the experimental θ = 2/3ML reconstruction does appear favourable.
A calculated DOS for each of these selected reconstructions is shown in figure 6.11. The DOS shown is in the surface-most Si and C atoms together with the potassium, i.e. the bulk DOS is not shown in these plots.
With higher temperature MD studies and simulated annealing, the overall system reached as low a minimum as the 300K anneal studies shown in figure 6.9. The implication is that the quench forced a reconstruction that is not entirely physical. An annealing study of the experimental surfaces showed significant reconstruction of most surfaces and, physically, more correct results. For accuracy, both are represented in figure 6.12, but the annealed results are what is reported.
Figure 6.13 shows the final lowest energy configurations found for the 1ML, 2/3ML, 3/4ML and 5/6ML coverages, respectively. As can be seen in the graphs in figures 6.9 and 6.10, the topologies found for these coverages were more stable than the previous interpretations of these surfaces, and maintain a similar electronic structure from figure 6.11.
The results of the calculations presented here show many interesting properties. Firstly, these results differ significantly from those of Derycke et al. for the θ = 1ML surface. The reconstruction found here to be most favourable is shown in figure 6.13(a). It was also concluded that the semiconductor to metallic behaviour detected previously by experiment STS appears to be entirely due to the potassium overlayer.
For lower coverages, the surface described experimentally for the 2/3ML coverage appears to be consistent with the computational results presented. However, the total energy of the system for this coverage is remarkably close to two or three other surface reconstructions that were simulated. MD of the 2/3ML surface described by SXRD showed the dimers twisting into a reconstruction more closely related to the final structure for the 1ML system proposed in this work.
Concerning the LEED results reported in the literature, the 1ML potassium positions proposed by Derycke et al. is based on a straightforward translation of the LEED into an atomic model. However, the results from calculations clearly show that this surface is energetically quite unfavourable. The lowest energy structure for the 1ML found by calculations is 2 × 4, (or c(2 × 4)). Disorder in the “× 4” direction might yield a 2 × 1 LEED, for example, due to the existence of different domains and corresponding domain walls in that direction. It is possible that the 1ML surface is a limiting case that is not fully realised in practice, because the potassium atoms are “too stressed”, that is, they are forced to be too close to each other and they relieve that stress by the formation of domain walls, so that coverage is actually slightly below 1ML under experimental conditions. That disorder in the “× 4” direction would yield a 2 × 1 LEED. Derycke et al.’s 2/3ML, which is very close energetically to this work’s 2/3ML surface topology, explains the experimental LEED results well. However, the 2/3ML coverage surface from computational results is an R30°(2 x 3), which is far more difficult to explain, based on experimental LEED results. Given the very small energy discrepancy between the 2/3ML from Derycke et al. and the computational surface, the 2/3ML topology proposed by Derycke et al., is viable, based on the computational results discussed in section 6.23.
The DOS results shown in figure 6.11 show some very interesting properties. The silicon carbide directly under the potassium surface is metallic. The semiconductor/metallic transition shown by Derycke et al. with increasing potassium coverage is based on STS results. This is explained quite well by the computational results showing the closing of the gap in the potassium DOS with coverage.
Calculations for this work have also shown that where potassium alone is allowed to coalesce, in the absence of a substrate, in the same arrangement as that found in the presence of SiC, the potassium can form a stable crystal structure. This is depicted in figure 6.14. Although the potassium forms a similar shape, the potassium-potassium distances are larger than in the presence of the substrate. However, when potassium is placed at a greater distance from the surface it does migrate towards it.
The implication is that the potassium is physisorbed onto the SiC surface, and that the bond between the two layers is nothing more than physisorption. This is also the interpretation put forward experimentally.22 An examination of the charge transfer between the potassium and the substrate sheds light on this matter more clearly, and a table of the charge migration is in Appendix III. The chemisorption energy was estimated at -1.05eV for the 1ML SXRD derived surface, -1.80eV for the 2/3ML coverage, whereas for the surface reconstructions found here, a chemisorption energy of -1.38eV and -1.83eV is found for the 1ML and 2/3ML coverages, respectively.
The results presented here agree with the experimental results as to what constitutes saturation on this surface, and in the fact that potassium atoms preferentially tend to reside in these so-called “pedestal” sites. When original atomic positions are placed off-pedestal, even in energy minimisation without MD reconstruction, the potassium atoms “fall into” these pedestal sites. However, after MD, this is not necessarily the case, especially for higher coverages. The final lowest-energy reconstructions are shown in figure 6.13. It was seen that once potassium coverage reaches about 5/6ML, the potassium atoms are no longer inclined to remain in the previously described “pedestal sites”. The experimental evidence supporting the pedestal site placement of potassium is not conclusive, based on ref []. The proximity of the potassium atoms to one another appears to force the potassium atoms out of these sites, but they cannot escape far enough from one another to prevent shells overlapping. Overlapping shells appear to be the mechanism by which metallic behaviour is then observed.
The final reconstructions shown in figure 6.13 show a remarkably different topology from experimental interpretations. At a coverage of 5/6ML, the potassium atoms cease to remain in the “pedestal” sites. The electronic structures of the 2/3 and the 1ML coverages appear to be similar to those reported in the literature.
Appendix III contains more detailed analysis of this system. This research is the subject of a pending paper submitted to Applied Physics Letters on the 25th of March 2011 for all coverages and comparison with the SXRD results of Derycke et al.. A paper specifically discussing the dimer reconstructions was submitted on March 25th 2011 for the “Sankey Festschrift” edition of Physica Status Solidi B, which is a special edition in tribute to Prof. Otto Sankey on the occasion of his 60th birthday. A presentation was delivered on this work at the APS March Meeting in Dallas in March 2011.
This chapter discussed the potassium-induced surface reconstructions on 3C-SiC(100). Owing to the large number of possible configurations of potassium on the surface, FIREBALL is very well suited to problems of this nature. Due to the fact that, within FIREBALL, all interactions are pre-calculated, it was possible to simulate over a hundred different surface reconstructions and evaluate their respective merits.
These results show a new model for the surface reconstruction. The reason for the semiconducting to metal transition on the surface was explained by computational techniques. Chapter 7 describes the study of another surface in silicon carbide and a semiconducting-insulator transition.
Following on from the metallisation of silicon carbide studies in Chapter 6, this chapter discusses the work carried out on the 3C-SiC(111) surface and the equivalent surface of the 6H polytype, the (0001) surface. In this chapter the metal-insulator Mott-Hubbard (MH) transition is discussed and an alternative mechanism to that reported to date in the literature for the transition is proposed. Section 7.1.1 discusses the 3C-SiC(111) and the 6H-SiC(0001) surfaces, for which the crystal types were introduced in section 6.1. Section 7.2 describes the experiments carried out. Sections 7.3 and 7.4 provide results and conclusions.
There are a large number of reconstructions presented by the (111) surface in 3C-SiC or its equivalent in the 6H polytype, the (0001) surface. Regardless of the polytype, the (111) or (0001) surfaces show the same reconstructions. The carbon-terminated surface has been shown to exhibit a (2 × 2) and a (6 × 6) reconstruction. The silicon terminated surfaces include the so-called silicon-rich surface (3 × 3) and the silicon-not-so-rich surface (√3 × √3). [59, 85, 86] The silicon-not-so-rich (111) surface of the cubic SiC and its counterpart on the 6H, the (0001) surface, has been reported throughout the literature to exhibit a Mott-Hubbard transition. Theoretically, the surface should have a single surface band halfway through the bandgap. However, experimental measurements show no such band by STM,[86, 87, 88, 89, 90, 91] and show two bands which are not accounted for theoretically in the LDA or the GGA.[92, 93] Theoretically, the expected single band in the gap is splitting into two bands, known as a Mott-Hubbard transition.
This 3C-SiC(111) √3 × √3 surface and the 6H-SiC(0001) √3 × √3 surface were explored in detail using FIREBALL to further understand this possible Mott-Hubbard Transition. Using STM[86, 87, 88, 89, 90, 91] the √3 × √3 surface for both the 3C and the 6H polytypes have been shown to have a fully occupied state about 1eV above the valence band minimum, resulting in a bandgap of about 2eV. [86, 87, 88, 89, 90, 91] Conversely, computational results in the LDA and the GGA indicate a conflicting result of a dangling bond-related half-filled band within the gap.[92, 93] This discrepancy between experimental and computational results has been previously attributed to a Mott-Hubbard transition.[10, 94, 67] However, by taking advantage of the unique efficiency of the FIREBALL package, results are shown that present compelling evidence of a far more elegant explanation of the mechanism for this apparent band splitting. FIREBALL is ideally suited for problems that require long-run MD, of the order of seconds of simulation time. This is vital when exploring phenomena like soft-phonon interactions, which may be the reason for the band split.
In section 7.1.2 of this chapter, the Mott-Hubbard transition is explained in more detail, in section 7.1.3 soft-phonon interactions are discussed. Sections 7.2 and 7.3 discuss the actual work carried out with FIREBALL on this interesting surface to explore the band splitting.
The Mott-Hubbard transition occurs when a surface band is observed to split into two bands. In the case of the 3C-SiC(111) and 6H-SiC(0001) √3 × √3 surfaces which this work focuses on, the dangling surface Si bond within the bandgap splits into two bands. This occurs when the ratio of parameters known as the hopping coefficient, t, and the correlation coefficient, U, reach a critical point. The parameter t is the energy overlap, in eV, between two wavefunctions, and U is the inter-atomic Coulomb interaction in eV. The relationship between these two parameters and the bands is depicted in figure 7.2.
In figure 7.2, E0 is the ground-state energy of the band and E0 + U is the energy of the band plus electron-electron repulsion. Within ab-initio simulations, as the overall energy of the system is reduced (i.e. the Hellman-Feynman forces are minimised), there is one band which can split when the ratio t/U drops below some threshold . Such a split is not accounted for within the LDA or the GGA.
The dependance on the ratio t/U makes these two parameters the fundamental variables to calculate. Mathematically, this behaviour can be described based on the Hubbard model, which is computationally expensive. A full mathematical description can be found in refs [[93, 95]] or in textbooks on the subject (for example, ref []). The model is specified by the one-electron dispersion relation of the dangling-bond band ϵ(k) = ϵ0 + t(k) and the intra-atomic Coulomb interaction, U, which is the energy required to remove an electron from one atom and add it to another. When U is greater than the width of the band, the quasiparticle energies are:
The double-degenerate band in the one-electron model is replaced by two narrower bands that are separated by the energy U in the Hubbard model. The strong electron correlation effects open an energy gap which is proportional to U between the two bands E1(k) and E2(k) that can be either completely empty or fully occupied.
In the √3 × √3 SiC surfaces the crystallographic surface reconstruction has been well established.[10, 67, 97, 98, 99] The non-metallic surface states seen by refs [[91, 67, 6]] have been attributed to a Mott-Hubbard metal-insulator transition with the single dangling bond per unit cell √3 × √3 band splits into a fully occupied lower band and an empty upper band.[10, 91, 67, 97, 98] Refs [[91, 94, 100]] report a dangling-bond bandwidth of ~0.45eV for this surface and experimental work using angle-resolved photoemission spectroscopy, STM, and computational techniques places the value of U as being ~2.2eV. By the Mott-Hubbard metal-insulator transition, the splitting leads to two surface bands of about 0.2eV width. Mott argued that the transition is sudden, occurring when N1/3aH ~ 0.0220, 21(N is the number of electrons), which is equivalent to when t/U = 0.218 for the √3 × √3 t has been calculated at 0.05 making this ratio 0.25.
The unique nature of a phenomenon known to cause a metal-insulator transition was studied from a differing standpoint in this work. The hypothesis is that an MH transition is not the reason for the surface behaviour. Rather, the existence of soft-phonon interactions on the surface is shown, and by carrying out careful electronic structure and energy calculations, show how these soft-phonons may be the reason for the experimentally-observed “splitting” of the single-band in the gap.
Another explanation for the surface behaviour discussed above is what is known as a soft-phonon interaction. Such interactions have been successfully argued as explanations for novel behaviour[59, 85] in a number of systems, such as in the much-debated In-Si(111) surface which exhibits “quantum wire” behaviour.
Within the soft-phonon model, the band-splitting surface behaviour is due to high frequency repetitive translations on the surface. Such oscillations may be well into the teraHertz regime, where they are outside of the sampling rate of current experimental techniques. For example, scanning tunneling microscopy (STM) and atomic force microscopy (AFM) sample in the megaHertz range. Other very important surface techniques, such as LEED, are not time-resolved.
Consider, for example, a see-saw like motion, one atom of a correlated pair rises as another falls and vice versa. Such a coupling would allow for a soft-phonon - a near-zero energy loss or gain, yet rapid change in atomic position. This is in contrast to a phonon, which is simply a vibration with appreciable energy change throughout the overall system of interest. The zero change in energy means that soft-phonons are difficult to detect, but due to this motion, can have profound effects on the electronic structure. The conclusion that the results of this thesis suggest that just such an interaction is occurring on the √3 × √3 surfaces studied in this work. With such see-saw motion, the frontier electrons in the dangling bonds will appear shifted. In the case of experimental measurements of the surface density of states (DOS) the net effect would be the observation of two bands, essentially the “sum” over time of a large number of fast moving measurements. In contrast, computational techniques are usually carried out on a system in which the atoms have been made come to rest- by reducing the Hellman-Feynman forces to a minimum, for example. This means that the atomic positions are “frozen” during the “measurement” of the surface DOS, leading to a single band that is in the middle of the upper and lower “virtual” bands from experiment.
In favour of the argument that soft-phonons are responsible for this behaviour is the temperature dependence of the √3 × √3 surface for a metal-insulator transition to occur. As temperature increases, so does the occurrence of the surface phonons leading to the metal-insulator transition. By carrying out long MD simulations it should be possible to observe soft-phonon interactions and by detailed analysis of the electronic structure of the dynamical surface an explanation of the observed band splitting can be shown based on the soft-phonon interactions.
The mechanism itself has been used to explain many similar surface transitions in the literature. For example, the highly controversial In-Si(111) surface, which is discussed later in this thesis, Sn/Ge(111)-√3 × √3 to 3 × 3,[100, 101, 102, 9] the In-Si(111)-4 × 1 to 8 × 2[103, 104] and the SiC(100) surface.
FIREBALL is ideally suited for studying soft-phonons as they present themselves in very-long-run MD simulations containing a larger number of atoms. Such studies would be prohibitively expensive with PW based tools. The computational efficiency of FIREBALL, however, allows us to study (relatively) large supercells over extended numbers of time steps (>105 @ 0.25fs).
In the case of the √3 × √3 surfaces that are discussed in this chapter, it is important to look at the electronic structure of the surface. With an adatom-above-a-trimer configuration such as the √3 × √3, there is a single half-filled dangling bond from each surface adatom. This dangling bond has been previously identified as the contributor to the Mott-Hubbard transition. Within a soft-phonon model the same bond is seen as two bands due to its stretching as it oscillates on the surface between two states at very high frequency.
To show the existence of such phonon interactions, a well prescribed low-temperature model of the surface is required. Due to the intensity of study of SiC there is much information available (e.g. ref []). The low-temperature model of the surface is depicted below, in figure 7.3. Extensive MD must be carried out to see how the surface develops with time, and thus identify the dynamical processes involved.
FIREBALL has previously enjoyed much success in problems of this nature, for example the reversible phase transition in ref []. The computational efficiency of the code allows for both the search for differing energy minima[9, 103] and to carry out long run MD[100, 102, 104] In this study an optimised simple sp3 basis was used. In order to optimize the simple basis set, calculations were
performed for bulk SiC, Si and C and the basis set that yielded lowest total energies and best equilibrium lattice parameters was selected this is the same basis set optimisation described in [].
In testing and experiment, results of a bulk modulus of 4.42Å and 4.45Å were compared, chosen based on the results of Trabada et al. in ref []. Stress or strain may adversely affect results and using two bulk moduli avoids such issues. By developing a basis and supercell for each of these bulk moduli, errors due to inadvertently causing stress or strain within the supercell can be avoided.
It has been shown that stress or strain caused by a slight change in bulk modulus may adversely affect computational results. In testing and experiment, results using bulk moduli of 4.42Å and 4.45Å were therefore compared. By developing a basis and supercell for each of these bulk moduli, effects due to inadvertently causing stress or strain within the supercell can be eliminated.
A number of calculations were carried out using this setup. Firstly, low temperature phases were identified by quenching. Secondly, extensive MD runs were carried out. Thirdly, hundreds of surface BS diagrams were generated for the same number of phonon-induced surface configurations. Finally an umbrella sampling and weighted histogram analysis method was used to identify the potential barrier between surface states, this method is explained later in this section. For each of the two polytypes studied, molecular dynamics simulations were carried out. Initial simulations were carried out to identify if there was a soft-phonon distortion effect actually occurring on this surface, as well tests to confirm the lowest energy configuration. Because there is only 1 silicon adatom in the √3 × √3 surface configuration this was carried out in a 2√3 × 2√3 supercell consisting of five layers of SiC with the single Si adatom per √3 × √3 unit. The bottom-most layer was fixed and hydrogen-terminated to simulate bulk. In order to simulate as many time steps as possible within resource limits, a 2√3 × √3 and √3 × 2√3 supercell was generated for each polytype and bulk modulus combination. A surface of this size contains two silicon adatoms and therefore a supercell of this size contains the minimum number of atoms to be simulated that still allows for soft-phonon effects to be observed. To avoid discrepancies, larger supercell calculations were carried out to check that errors do not occur in the smaller cells.
Once a “steady-state” of oscillations was observed in each of the MD simulations, one full oscillation was selected from the atomic position files. For each time step in the selection, the band structure was calculated. By “summing” over the band structures, it was possible to show that the degenerate “split” bands emerge.
Finally, the UMBRELLA sampling method was used to get a positional density of states in real space for each case. The UMBRELLA sampling method is an application of the Weighted Histogram Analysis Method to the energy and forces of various configurations of the system, as it processes through an MD run. Simply put, it calculates the probability distribution of the configuration of the system. It is an excellent analytical tool for studying the most favourable positional configurations of the system in study.
The UMBRELLA energy and forces are calculated relative to the distance between 2 atoms. The form of the umbrella potential is the harmonic relation:
where K is the force constant, dab is the distance between atoms a and b, and d0ab is the fixed value given as the centre of the window.
The implementation of this within FIREBALL is based on the work of ref [] the reaction co-ordinates chosen for the UMBRELLA sampling in each case were the distances between one of the topmost silicon atoms and whichever atom is directly beneath it. In the case of the 3C, this is a silicon atom in the layer directly beneath the surface adatom of interest. In the case of the 6H-SiC, this is a silicon atom in the 4th layer. The sampling bins are set at 0.1 the total maximum displacement between the atoms at the chosen reaction coordinates.
The UMBRELLA sampled data, which is a series of parabolas representing the probability distribution of the sampled data within each bin is then passed into the WHAM analysis. An example of the UMBRELLA sampled data is in figure 7.11, the procedure is outlined by Roux. The UMBRELLA data is divided into smaller bins and the probability of these is used to measure the free energy in the system with respect to the displacement between the centers. The output graph gives us a feel for the local energy wells within the dynamic system. Graphing the free energy versus the displacement of the two centers shows where these potential wells are, and the energy required to pass between them.
The reaction coordinates chosen for both the 6H and the 3C-SiC were the surface adatom and the atom that is directly below it in the lattice. In the case of the 3C, this is in the layer below, however in the 6H it is three bilayers below. A representation of this is shown in figure 7.4, below.
For each of the 4 sets of experiments, namely 3C and 6H with a bulk modulus of either 4.42Å or 4.45Å an extended MD was run. Due to the large number of timesteps of these studies, only an excerpt of the z position of the surface adatoms is shown in figures 7.5a-f. In total sixteen long-run MD simulations were carried out for this study. The simulations carried out break down as follows- for each of the 4 categories, namely 3C and 6H-SiC with a bulk modulus of either 4.42Å or 4.45Å, there is a 2(√3 × √3) surface simulated, a 2√3 × √3 and a √3 × 2√3 and a √3 × √3. The √3 × √3 surfaces, which contains one surface adatom per study, are not applicable in the soft-phonon results here but are very useful when discussing the electronic properties later.
The 2√3 × √3 and the √3 × 2√3 surfaces were simulated because these two surfaces allowed us to explore far more MD steps than the larger surfaces due to the lower number of atoms in the supercell. For these surfaces, four k-points were chosen for use during MD of the 2√3 × √3 and the √3 × 2√3, two for the 2 ×(√3 × √3) and eight k-points for the MD simulations of the √3 × √3. The choice of k-points meant that calculated energies could be directly compared between supercells of different surface sizes.
Figure 7.5 shows examples of the soft-phonon motion, two coupled adatoms distance above the surface are shown plotted against time, for these graphs where each time step is 0.25fs. The correlation between the coupled adatoms is clearly shown in these pictures. The period of these oscillations is typically about 1000 time steps, which corresponds to a time of 250fs or a frequency of 4THz, well outside resolution of current experimental techniques. All graphs in figure 7.5 are from simulations at 300K.
The graphs in figure 7.5 were produced for less time steps over a range of temperatures from 50K to 1000K in the MD. The period of oscillations reduced with reduction in temperature, down to 400fs for the 50K studies and up to 200fs for the 1000K studies. These periods correspond to frequencies of 2.5THz to 5THz.
A “check” set of experiments was carried out in all cases graphed in figure 7.5. The check was that one of the surface adatoms was fixed in position in the 2√3 × √3 and the √3 × 2√3 studies, or two surface adatoms in the 2×(√3 × √3) studies. The free surface atoms oscillations were severely reduced in amplitude by comparison to the above graphs. The fixing of one atom out of a pair that oscillate in tandem to one another allowed a confirmation that there is a dependance of one adatom in the pair on another.
Observing the existence, frequency and coupling of these soft-phonons is not enough to explain the surface behaviour. These results confirm the existence of the proposed soft-phonon and a closer look is now required at the electronic states of the systems. In order to do this, a band structure was generated for a large number of timesteps for the studies graphed above. One hundred surface bandstructures were calculated for every 300K simulation within one period of oscillations as seen in figure 7.5. A summed bandstructure for the 2×(√3 × √3) studies is shown in figure 7.7.
In order to carry out the band structure calculations, 41 “special” k-points were chosen along the axes depicted below (figure 7.6) where the vertices represent the high symmetry points in the 1st Brillouin Zone. The band structures shown in figure 7.7 are for the 2 × (√3 × √3) surface, which has four surface adatoms and therefore four surface dangling bonds. The oscillatory nature of the surface can be observed by the spread of the bands in the gap about the static band positions. Further calculated bandstructures are in Appendix IV.
Calculated Local Spin Density Approximation results for these two surfaces are shown in figure 7.8, this figure also shows LDA results for the energy-minimised surface for comparison. As can be seen, the observed summed-bands in figure 7.7 closely approximate the LDSA and experimental results shown in figure 7.8.
In order to depict the rapidly changing band structure with respect to a dynamical surface, figures 7.9 and 7.10 show a three part description of the soft-phonon oscillation. The first part is a the z position of each surface atom with time. Directly below is a graph of the total energy of the system over the same time and, thirdly, the band structures are shown for specific points in the time-dependent graphs. For this study a single oscillation from the Z-data graphs was chosen, i.e. an approximate run of each atom pair going from low to high and back again. This single oscillation is then graphed in terms of total system energy versus adatom height in the z direction in the graph. The band structure at the peaks and troughs in this graph are then shown.
The final set of studies of the surfaces was carried out by applying the UMBRELLA method. As described in Section 7.2, choosing a spring constant that results in well defined and meaningful curves is of utmost importance. For each of the datasets, some few hundred spring constants were tested for each of the studies and WHAM was applied to those that resulted in well-defined curves. An example of the UMBRELLA output is depicted in figure 7.11. The resultant WHAM graphs are discussed after a description of the UMBRELLA results.
A k constant of 35 was decided upon for the 6H-SiC and a k constant of 120 was used for the 3C study, chosen by trial and error. This resulted in the graph in figure 7.11 and similar graphs for the other supercells. Applying the UMBRELLA results, of which an example is shown in figure 7.11, to the WHAM routines results in a curve showing the free energy per mole versus the distance between the coordination numbers. This measurement of the free energy is only qualitative, however, as can be seen in figures 7.12 and 7.13, there is a potential barrier which is also quite small. Although graphs 7.12 and 7.13 are not quantitative, it is clear that the potential barrier outside of the ideal windows goes to infinity on either side of the small hump. This implies that there are two metastable surface states between these coordination numbers, indicated by the wells on either side of the small barrier and that outside of the region the potential tends to infinity.
These results indicate the presence of soft-phonon motion being the likely cause of the metal-insulator transition observed on this surface. Once the system is in motion, it visits metallic and semiconducting states at high frequency, causing the observed two-band formation that has been experimentally observed.
By comparison to other cases where soft-phonon interactions have been observed [10, 76], there are many similar features in the behaviour of the present systems. However, unlike other systems, there has not been any evidence in the literature for a change in the surface reconstruction on the SiC(111) and SiC(0001) surfaces studied here. The observed oscillations, being confined to one dimension in general, make it difficult to observe a change between metastable states in experimental conditions and this may be the reason that no surface modification has been discussed in other works. The results of “freezing” one or more surface atoms and allowing the others to move freely provides further compelling evidence that the position of the atoms is dependent on the frontier orbitals (the dangling bonds) as the unfixed atom moves far less when this is done.
Compelling evidence for the soft-phonon hypothesis can be seen in figures (7.11-12). The energy and band structures at the lower energy points are sufficiently varied to show the band splitting that was previously supposed to be a Mott-Hubbard transition. The WHAM results further support this, showing a pair energy minimum and a very low potential barrier between the two metastable states.
Figures (7.4 a-f) show the frequency of these oscillating states- in the teraHertz range for both the 3C(111) and the 6H(0001) surfaces studied. This is well out of the range of current experimental techniques. The range of the oscillation amplitude is also very small, in the region of hundredths of angstroms. The effect on the band structure is profound, as shown in figures 7.9 and 7.10.
Mott-Hubbard transitions support much of the properties observed on this surface. However, such a transition is often used as a catch-all term for certain phenomena, while this present study shows the possibility of a far more elegant explanation of the systems’ behaviour. This adds up to a few points. There are two energetically degenerate ground states of these systems. At room temperatures, the system dynamically jumps between the ground states. This is outside of the time resolution of modern experimental techniques. The only current way to observe this is via MD simulations.
This chapter describes the surface transition experimentally observed on the (111) surface of 3C-SiC and the (0001) surface of the 6H-SiC. This transition has previously been attributed to a Mott-Hubbard transition, however it is argued that the same transition can be described by a soft-phonon interaction. It is shown that the energy required to oscillate between two surface states is remarkably low and as a result a soft-phonon may well be the mechanism for the transition.
In the next chapter, future work is discussed on another surface that shows a phase transition, the In-Si(111) surface. This surface has been shown, computationally, to exhibit a soft-phonon transition. There is also much experimental data available on this surface, including many RA spectra. With careful development a calculated RA spectrum based on a dynamical surface would offer strong evidence in favour of the soft-phonon theory over other prevailing ideas about the In-Si(111) surface.
In Chapters 6 and 7, two separate phenomena of nanoscale systems were discussed, namely semiconductor to metal transitions and the Mott-Hubbard transition. In this Chapter another type of behaviour is discussed- nanoscale “quantum wire” behaviour. This occurs when a metal overlayer on a substrate forms chains that are far longer than they are deep or wide forming a quasi-1D system. This chapter discusses what may be a dynamical surface, the indium induced 4 × 1 reconstruction on Si(111) surface, which has been described as exhibiting quantum wire behaviour. As discussed in section 8.1.2 of this chapter, the behaviour may be due to soft phonon effects. Preliminary computational work has been carried out on this surface, as described in section 8.3. In section 8.4 an experimental surface analytical technique called reflectance anisotropy spectroscopy (RAS), which has been important in the study of this surface, is discussed. Section 8.5 describes how this technique can be modeled in order to compare to experimental results. Section 8.6 discusses the work carried out so far in pursuit of the goal of modeling RAS for a dynamical surface.
So-called “quantum wires” are a quasi-1D construction atop a substrate, with the electrons confined in the other two dimensions. The term wire is derived from the fact that the conductor has length but negligible depth or width. There is much technological interest in this phenomenon, as it could lead to nanoscale wiring systems of high conductivity.
As a result of the the electron wavefunctions confinement in two spacial dimensions, unusual properties present themselves and physical restriction or quantum confinement can occur. These properties may manifest themselves in unusual effects that depend on electron behaviour such as the optical and magnetic properties of the systems in particular.
The novel properties of these systems may allow us to engineer devices that are applicable to fields such as nanoelectronics, spintronics or other applications. By understanding the mechanisms by which certain unusual properties of these systems arise, it may become possible to tailor-make systems with near-ideal properties to a desired function.
Prototypical of such quantum wire behaviour are variants of the In-Si(111) reconstruction. This surfaces has been discussed in the literature as a subject of interest since the 1960’s. The known reconstructions are currently of intense research interest and discussed later in this chapter.
As part of the studies for this project, it was planned to carry out an ab-initio RAS calculation of this surface. However, due to restrictions with respect to the rewriting of the code, this is now in “future work”. Here, preliminary results on the In-Si(111) surface are discussed, some of which were presented to the American Physical Society at the March Meeting in 2008 in New Orleans, have also been discussed in poster presentation form at the Institute of Physics Spring Weekend in Birr Castle, Ireland in March 2006 and at the 2nd National Conference on Condensed Matter & Material Physics conference in April 2007 in the University of Leicestershire, U.K.
RAS is an excellent experimental technique for analysing quantum wire systems. It is a surface analytical technique that measures the anisotropy in the near-normal incidence reflectance of linearly polarised light along two orthogonal axes in a material surface. These data are highly surface-specific and have been shown recently to provide an optical signature of conducting nanostructures.
In this chapter, section 8.2 discusses the In-Si(111) system and computational findings carried out for this work and results are compared to similar computational and experimental results. RAS as an experimental technique is discussed in more detail in section 8.3 and the way in which a routine RAS calculation engine could be developed for the FIREBALL suite is discussed in section 8.4. Section 8.5 discusses the current progress to this end in FIREBALL. The In-Si(111) system was chosen as an ideal starting point due to its predominance in the literature as it has been the subject of intense study in recent years.
A large number of reconstructions has been observed on the In-Si(111) surface with coverages of 0.5 monolayers. The first details about this surface in the literature are from 1964 in work by Lander and Morrison.1 In the last decade, a lot of work has been published on the surface reconstructions presented by this system following Yoem et al.’s 1999 and 2002 work which showed that this surface not only exhibits one dimensional metallic properties, but also a temperature induced reversible phase transition at 120K, leading from the 4 × 1 to an 8 × 2 reconstruction. Initially, for the 4 × 1, various models have been proposed by Sarinin et al.[111, 112] and Zotov et al., Kumph et al., Bunk et al.,[115, 116] on the basis of results from various techniques such as RAS,[117, 5] STM,[112, 5, 118] LEED, PES, XRD and HREELS.
In 1997, Kraft et al. studied the various reconstructions using STM and STS, identifying numerous “striped” structures on the surface at sub-monolayer coverages. This work resulted in the phase diagram in figure 8.1. Following that paper much interest in the so-called “striped structure” and the 4 × 1 began.
The room temperature phase of the In-Si(111) surface exhibits a 4 × 1 reconstruction, the generally-accepted topology of this surface is the Bunk model, which came about after a period of deliberation between various proposed models.[111, 112, 113, 114, 115, 116]
Figure 8.2 shows the bunk model surface reconstruction. It consists of two indium chains separated by a silicon chain. The reconstruction itself is similar to the Pandey type reconstruction for Si(111)-2 × 1 surface. One interesting result of this effort in the In-Si(111) 4 × 1 is that, due to a large energy gap between the wire and the substrate, it has been shown that the indium chains appear to carry charge without any coupling or dependency on the surface, and thus may be the smallest known wires possible. The 4 × 1 surface itself appears to be a 1-D conductor, with electron confinement in the other two dimensions and is a Peierls-type conductor.[109, 110, 124, 125] In a Peierls-type conductor the electrons do not behave as they do in a Fermi liquid, instead they exhibit “Luttinger liquid” behaviour, which is described in detail in reference [].
Published computational work on this surface agree with the Bunk model for the 4 × 1 reconstruction. The first of these theoretical works to appear in the literature was by Cho et al., where the conclusion was drawn that the Bunk model for the 4 × 1 was confirmed to be the lowest energy reconstruction. The conclusion also supported experimental data from angle-resolved photoemission, scanning tunneling microscopy images and core-level photoemission data. Other computational studies published include Nakamura et al., and R. H. Miwa and G. P. Srivastava, who independently confirmed the Bunk Model via PW in both the LDA and the GGA and tight-binding methods in both the LDA and GGA. These studies agreed with Bunk’s 4 × 1 model as being the lowest energy reconstruction at this temperature by total energy calculations from minimisation of the Hellman-Feynmann forces and evaluating the overall structural change. Around this time Tsay et al. computationally evaluated the three phases of the surface in the LCAO using LDA and evaluated the total energies. 
The LT reconstructions are still hotly debated in the literature. In this regime, there are two identified phases, a 4 × ”2” and an 8 × 2. The phase transition 4 × 1 -> 4 × ”2” -> 8 × 2 has been attributed to a CDW driven by a single band[109, 110] or triple band Peierls instability, or that it is due to a periodic lattice distortion,9, 23 or a many-body interaction effect, or some combination of these. The idea of the phase change being due to a CDW driven by a Peierls instability was questioned by Kumpf et al. based on their SXRD results and they attributed a 3D CDW. The Peierls deformation is driven by a charge density wave (CDW) as described in Peierls book on the subject.
Furthermore, there has been no general agreement on the structure of these LT surfaces, with models proposed suggesting the 8×2 surface is metallic, but with a lower density of states at the Fermi level,[110, 132, 133] semimetallic, and semiconducting with a fundamental energy gap of 0.1–0.3 eV.[134, 135, 136] Cho et al. have concluded that the phase transition was not due to a CDW. However their interpretation of results is not completely accepted. Tsay proposed a candidate reconstruction for the 8 × 2 surface in 200527 which corresponded with Kumph et al.’s findings. This is shown in figure 8.3, this is the “trimer model” for this surface.
The 4 × ”2” reconstruction remains elusive. As figure 8.4 indicates, LEED results do not show sharply-defined spots indicative of perfect order in the x 2 direction, but instead show “streaks”. Again, this has been the subject of many ongoing debates. It is difficult to accept, based on this “streaked” LEED pattern, that a CDW would be the driving mechanism for the phase-change as a CDW is a theoretically macroscopic feature across the whole surface. Therefore a CDW-driven shift would affect the entire surface simultaneously, and not give rise to the non-uniformity seen in the LEED for an appreciable amount of time.
The effort to understand the phase change and streaked LEED pattern has been considerable. With the possibility of quantum wire behaviour in the 4 × 1, the subject of the LT phases is of significant technological importance.
Two excellent methods for probing such anisotropic electrical properties is the use of polarised optical techniques such as RAS and Raman spectroscopy. The first RA spectra of the In-Si(111) reconstruction were by Pedreschi et al. in 1998. Later, Fleischer et al. and Chandola et al. succeeded in achieving a RA spectrum for lower energy incident light and over the phase transition itself. These studies which reported on both RAS and Raman spectroscopy results, indicated both one dimensionality and that through the 4 × 1 -> 8 × 2 transition there is a split in the anisotropy at about 2eV, as well as 1.4 eV and 0.7eV. By using vicinal silicon and the technique of J. Viernow et al., all three succeeded in growing large domains of these 1D structures and then probing them experimentally. This then simplified the process of cooling the substrate and looking at the change in the RA spectrum as the surface altered from the 4 × 1 -> 4 × ”2” -> 8 × 1. The result of these evolution RA spectra was that subtracting one from the other gave further insight into the phase change. Recently, Chandola et al. concluded that a hexagon model for the LT surface was more likely over Fleischer et al.’s trimer model.
For the 4 × 1 structure, there have been a few successful theoretical RAS spectra reported,[126, 104] however theoretical RAS papers have yielded no further results for the LT phases. Due to the inherent limitation with DFT that the bandgap is underestimated, computational RAS results are relatively qualitative in nature. However Wang et al. succeeded in theoretically reproducing the overall features of the 8 × 2 and the 4 × 1. Although Wang’s results discussed the 4 × 1 Bunk model and their support for it from their results, conclusions about the 8 × 2 or the 4 × 2 topology were not forthcoming.
Using GGA with a PW basis set, Cho et al.’s analysis of the atomic structure of the SXRD-derived trimer structure yielded the result that the optimized 8 × 2 or 4 × 2 reconstructions were only ~8 meV/(4 × 1-unit cell) more stable than the 4 × 1 structure. The 4 × 1 structure that they found was in good agreement with ARPES experimental results, adding confidence to the 8 × 2 results, however 8meV is a very small difference and well within the error margins of these computational methods. In contrast with experimental results, they also found the 4 × ”2” to be metallic. This indicated that the processes’ responsible for the phase transition on the surface may have be something not included in their DFT calculations for the 4 × 1 → 4 × ”2” transition. Based on these results, Ahn et al. suggested with results from angle resolved photoemission studies that the phase change is due to a triple-band Peierls instability. Park et al.  suggested also that the phase change may have something to do with defects, and using STM observed the charge ordering to show the CDW nature of the ground state. Based on the interpretation of their STM and DFT results, Lee et al. argued that the LT 4 × “2” phase is not a band insulator, but is stabilised by many-body interactions. In 2005, Gonzales et al. showed that soft sheer motion may be responsible for the metallic behaviour and the phase change. Gonzales showed that there are 4 hexagonal 8 × 2 reconstructions into which the LT phase can form, rather than Kumpf’s trimer model. In higher temperature regimes these are shown to average out into the currently known 4 × 1 reconstruction by rapidly changing sheer motion interactions. Gonzales argues that this model (known as the GFO model) answers the issue of the 4 × ”2” streaks as well as describing the reasons for its metallic behaviour at RT. The GFO model of the system essentially states that there are four differing ground-state LT phases of the In-Si(111), each with 8 × 2 symmetry. Once the system reaches the critical temperature of ~130K all four exist in a constant state of flux on the surface due to soft phonon interactions, leading to an experimentally resolved 4 × 1 surface. In principle, this is a similar macroscopic behaviour to the small microscopic behaviour argued in Chapter 7 regarding the 3C-SiC(100) surface to explain the Mott-Hubbard transition.
In a letter, Yoem aggressively disagreed with Gonzales et al.’s conclusions, citing that dynamical fluctuations of this nature would be detected in photoemission experiments. Gonzales et al. responded quickly with another latter, showing that the model they propose is consistent with peak broadening as observed by experimentalists and that their model is consistent with all other available data. They showed that a soft-shear motion is not solely responsible for the reversible phase transition (RPT), but is definitely a factor in it. Very recently Gonzales et al., went on to explain the RPT in terms of a two-band Peierls transition and a soft sheer motion action rather than the previously thought 3 band Peierls transition.[141, 142]
Fleischer et al., succeeded in 2005 in acquiring an RA spectrum for lower energy incident light than previously possible. Their study indicated that soft shear phonons may well be the underlying mechanism. Recent work by Chandola et al., shows that the RAS and surface ellipsometry in the mid infra red can support one, but not all, of the energetically degenerate structures proposed by Gonzales et al.. However, using the PW approximation Cho et al. had previously shown that the hexagon structure is unstable. Yet recent positron diffraction studies support the hexagon model for the LT phase. In further contrast to these results, Wipperman et al. reports that the Gonzales model supports calculated results for quantum transport whereas any trimer model does not.
Very recently, Gonzales et al. has shown that theoretical STM results based on the dynamical model of 4 × 2 hexagons forming the observed 8 × 2 and 4 × 1 reconstructions match well with experimental STM results. This is further supported by Chandola et al.’s recent paper using a PW code that shows that experimental and theoretical results for a hexagon-type model are consistent in RAS results for a small supercell and for at least one of the possible 4 × 2 reconstructions.
It was proposed that by using the massively scalable computational abilities of FIREBALL coupled with theoretical RAS or similar ellipsometry results, it may be possible to further develop the model and explain some of the newer experimental results based on the dynamical surface model postulated by Gonzales. That is, to calculate a RAS for all the Gonzales structures and see if their interaction can show a minimisation of the 1.4eV peak when summed together.
In terms of actual computational work carried out on this surface, an Fdata and molecule file have been developed. Gonzales et al.’s soft phonon interaction required MD simulations, however, a quench has been carried out. The Fdata consisted of an LDA basis set with silicon and indium. The cutoffs for the silicon were 4.8Å for the s shell and 5.3Å for the p shell. For indium the cutoffs were 5.4Å and 6.3Å for the s and p shells, respectively. The basis set cut offs were default values for silicon and hydrogen in this case. The molecule file contained 20 indium atoms atop a slab of 350 silicon atoms, this represented a substrate of 9 layers of silicon and a surface size of 2×(4 × 1) in the “×1” direction. Figure 8.5 shows the supercell used for this calculation.
An energy minimization, surface DOS and bulk DOS has been generated for this system.
In the quench, the supercell reconstructed minimally. An average atomic movement in the indium of less than 0.1Å was observed. This implies that the structure is, initially, in a relaxed state. A table of charge migration and the DOS in the surface and in the bulk is discussed in this section.
The indium atoms closest to the surface silicon chain lost, on average 0.01eV, whereas indium atoms 11 - 20 gained a charge of 0.04eV. The charge lost on indium 1-10 was somewhat gained by the surface silicons. Also, silicon atoms 31 - 35 (not pictured), which are directly below the indium atoms numbered 11 - 15 had a similar charge transfer result to silicons 21 - 30. The other silicon atoms on the layer sharing silicons 31 to 35 are bound to the surface silicon atoms (numbered 21 - 30). Below that the change in charge on silicon atoms was of the order of 0.001eV. This implies some bonding is occurring between the indium and surface-most silicon atoms.
Figure 8.7 shows the calculated surface and bulk band structures. These were calculated with an 8 k-point Monkhorst-Pack mesh. The Fermi level for the slab was found to be -2.83eV.
As mentioned in section 8.2, these results are preliminary, however it is clear that there is some charge interaction between the indium and surface silicon atoms. The DOS in the bulk shows semiconducting behaviour, as expected and the surface DOS is metallic.
These results are an ideal starting point for the development of a RAS calculation from within the FIREBALL package. The rest of this chapter deals with RAS itself, in section 8.4 and how this can be built into FIREBALL-lightning as in section 8.5. Section 8.6 finishes off with the implementation of code to calculate optical absorption in FIREBALL as another stepping stone towards a complete ab-initio RAS calculation.
RAS is an optical technique for probing a surface. It is non destructive, and is relatively inexpensive to implement experimentally. As a result, since its inception in the 1980s, it has become an important tool in analysing surface reconstructions. The importance of RAS in that it excels as a technique when probing surface reconstructions that may have structures aligned in one direction, especially on a substrate that has cubic symmetry.
Typically, a RAS study is combined with one or more of a range of surface techniques such as scanning tunneling microscopy (STM), low energy electron diffraction (LEED), x-ray photoelectron spectroscopy (XPS), ultraviolet photoemission (PE) and inverse photoemission (IPE) spectroscopy. Combining these techniques with RAS provides insight into the atomic and electronic structure and the morphology of surfaces.
The principle of operation of RAS is that it measures the difference in reflectivity of normally-incident plane polarised light between two orthogonal directions in the surface, normalised to the mean reflectance:
where the reflectances r are complex Fresnel reflectance amplitudes. Most optical probing techniques are not very surface specific, because they penetrate many layers of the sample. However in the case of RAS, assuming the substrate is some bulk crystal with cubic symmetry, then the result is only surface contributions and negates the response due to the substrate.
Although RAS is applicable in any environment where the surrounding medium is optically transparent, it is usually used in ultra high vacuum (UHV) due to considerations for the specimen. In the last decade or so, RAS has been established as a powerful probe sensitive to extraordinarily low surface coverages such as ~0.01ML. Outside of its established use as a semiconductor and metallic surface analysis tool, RAS has been applied to remote surface stress/strain measurements,47 monitoring LCD fabrication and characterisation of thin films.
Experimentally, RAS may be regarded as a development of spectroscopic ellipsometry(SE),[150, 151, 152] a technique that can be used to evaluate the dielectric function of materials without the complexity of Kramers-Kronig analysis.[153, 56] The essential difference with SE is that linearly polarised light is incident on the sample close to the Brewster angle, whereas in RAS uses near-normal incidence of light.
In application, RAS requires an understanding of the Jones matrix and its applications[e.g. Hecht] in order to correctly set up the apparatus [See, e.g. , and ].This allows an experimental physicist to mathematically evaluate the effect the apparatus components has on the overall result. There are many examples of excellent discussions of this and the practical setup of the experimental apparatus in the literature [e.g. [153, 107] and references therein]. The essential use of the Jones matrix in RAS is that it describes what happens during an experiment.
There is no generally described method for interpretation of an RA spectra. This is especially difficult with single crystal surfaces since the response of the surface under investigation depends upon the complex dielectric function of both the surface region and the bulk. However, a useful method for analysing a RAS result is the “three phase model” developed from Fresnel theory by McIntyre and Aspnes, with extensions by Selci et al., Cole et al. and Del Sole.
A surface that has reconstructed may be considered to be an interface between two crystals with dielectric tensor components ϵixx, ϵiyy, and ϵizz. ϵixx is the x response of the dielectric tensor to x components of the incident beam in medium i, where i = 1(surface or thin film) or 2(substrate or bulk material). If the bulk material or substrate is isotropic, then the tensor components are the same along each material axis. This allows ϵ2jj, where j = x, y or z, to be replaced by the single component ϵb, the isotropic bulk dielectric constant.
The three phase model is the division of the surface into three distinct regions- bulk, thin film and vacuum. Using this model and the Fresnel equations, expressions may be derived for the s and p polarised light reflected from a surface. As RAS is only concerned with the reflectance of p-polarised light resulting from a p-polarised incident beam, only that equation need be quoted for experimental analysis of RAS.
In the normal incidence case, θ0 = 0 and this equation simplifies to:
where θ is the angle of incidence, and ϕ is the angle of polarisation with respect to the normal of the surface plane.
then substitution into the rpp (θ, ϕ) equation above yields the following expression:
where d is the depth of the surface reconstruction (medium 1) and λ the incident light wavelength and d must be << λ. This equation can aid in the calculation of surface dielectric tensor components and a general understanding of the RA response. This summarises how experimentally derived RA spectra can be analysed, as it converts the experimental quantity Δr / r into a material property (ϵxx - ϵyy = Δϵsurface). There is more information that can be gleaned from such procedures, and these can be found in any of the references in this section. (e.g. Weightman.).
In section 8.2, the advantage of calculations for comparison with experimental results when exploring the In-Si(111) surface. Calculating an RA spectra from first principles has allowed for the direct comparison of experimental models for surface reconstructions as well as giving insight into underlying mechanisms for behaviour in a system. The dependance of RAS on surface geometry means that a sufficiently accurate simulation technique to produce theoretical RA spectra from first principles is a highly desirable tool when for exploring surface topology models.
Building on the three phase model discussed in section 8.3 and the comprehensive earlier work put forward on matter-light interactions,[155, 163] Del Sole et al. developed in 1981 a generalised expression for the surface contribution to the reflectance,
where ϵb is the bulk dielectric function, and
where ϵij (ω; z, z’) is the non-local macroscopic dielectric tensor at the solid vacuum interface as described in depth in reference .
Early attempts to generate a computationally derived RA spectra were grounded in empirical methods. By fitting the the dipole tensors within the surface to simple harmonic motion, a qualitative picture is drawn up. A method like this was carried out by Wassermeier et al. in 1991 which resulted in qualitative results that were a good match with experimental findings.
In 1990 Manghi and Del Sole showed that a simple expression for the surface contribution to the reflectivity can be derived from eq. 8.12 above provided that two criteria are met. Firstly, that a supercell is sufficiently large that the surface and surface-modified bulk wavefunctions are described and secondly, that the off-diagonal terms of the dielectric tensor are small compared to the diagonal ones. The result is
where, αiihs with i = x, y is the dielectric tensor component of averaged half-slab polarisability. In principle, eq. 8.13 contains all surface contributions to the optical reflectance. The actual calculation usually requires some further approximations. The effort to solve αii by ab-initio methods is highly involved. The first ab-initio semiconductor RAS calculation was published by Morris et al. in the mid nineties within the LDA using a scissors operator to deal with the band gap issue.
Adolphi et al.  and references therein describe a method within the LDA for calculation of the αii term. The method described is still used in most calculations (for example in references [170, 171, 172]), however the application of TD-DFT and better XC approximations have been shown to yield a modest improvement in results. The issue is that this improvement has also been shown to be computationally far more expensive.
From DFT calculations in the LDA, the Bloch band eigenfunctions are obtained, with energy eigenvalues ϵn (k) in the first Brillouin zone, where n is the band index and k is the wave vector in the BZ. For a semiconductor, the Bloch states have an occupancy of 1 or 0, for valence band and conduction band, respectively. From Adolph et al., and references,[169, 174, 175, 176] the microscopic dielectric tensor in the optical limit of vanishing wave vector q → 0 is:
where v and c correspond to the case when n is in the valence or conduction bands. Ω is the volume of the reciprocal cell. η is a line broadening term.
Calculating αii in this way over a very fine real-space mesh to avoid numerical noise for the integration has been used in other implementations to yield an ab-initio RA spectrum. Results of integration algorithms tested within the FIREBALL package can be found in the previous chapter on the Create code.
This method has been directly used to determine the geometry of the 3C- SiC (001) 3 × 2 reconstruction and on the In-Si(111) 4 × 1 surface in previous studies. The ability to carry out such calculations routinely in FIREBALL would not only be a boon to its use as a tool, but would also allow for far more in depth analysis of the dynamical phases of the In-Si(111) surface discussed above. This in turn would not only yield a better understanding of this controversial surface, but also a greater knowledge of the underlying mechanisms of RAS, which are still debated.[178, 179, 180, 181, 182, 183, 184]
As a “first step” towards implementing the calculation of an RA spectrum, a new module was written for the FIREBALL code that would allow the calculation of absorption spectra. This new module was the first attempt to build a “plug in” for the FIREBALL-lightning code.
There are a large number of papers on the subject of simulation of an optical absorption via empirical or ab-initio methods. Within the LDA, the method used by Makowska-Janusik et al. was implemented. This has been applied to empirical calculations previously and is well described in the literature.
In order to accurately simulate optical properties of a system by ab-initio methods, a far more rigorous method should rightly be employed. In other words in calculating an optical absorption spectrum, the method outlined here is only weak approximation- there was no accounting for the electron-hole interactions, as would be done by the independent quasiparticle method, for example. The objective was an absorption spectrum within the LDA using the independent particle approximation that is at the core of FIREBALL. This project allowed for an evaluation of the philosophy behind FIREBALL-lightning and was not expected to result in a highly accurate outcome.
From Makowska-Janusik et al., the intensity of absorption can be expressed by the formula:
where I(ω) is the intensity of absorption of energy of frequency ω, j is the HOMO state in the molecule, k sums over the excited states and Γ is a Gaussian line shape widening parameter. Results from this module are not yet compiled.
Implementing equation 8.15 to calculate αii in FIREBALL-Lightning has a specific limitation. This issue, which was highlighted late in this work, is that the wavefuctions required in the equation require pure Kohn-Sham type wavefunctions. Kohn Sham is being implemented in FIREBALL, but had not been implemented at time of writing.
This chapter discusses a surface that shares properties with both of the surfaces discussed in Chapters 6 and 7, namely a metallic overlayer on a semiconductor and a dynamical surface. The natural progression of the work described in this thesis lends itself to the development of a method for calculating an RA spectrum for many different surface reconstructions within the soft phonon model for the In-Si(111) surface for comparison to experimental results.
In the studies on surface metallisation of silicon carbide, a new reconstruction for the 1ML was suggested. The surface model suggested fits experimental results well and has a total energy that is far lower than the SXRD-derived surface previously proposed by Derycke et al. The 2 × 1 LEED pattern reported from experimental results for this surface cannot be entirely explained by the computational surface model shown in figure 6.13a. It may well be the case that the overall pattern emerges as a result of a distortion in the “×4” direction of the c(2 × 4) underlayer, leading to the observed 2 × 1 due to potassium-induced boundaries interfering with the surface.
The computational results also suggest that the potassium is semiconducting up to and beyond the 1ML coverage due to its forced positions on the potassium surface, and that the silicon carbide surface is metallic throughout. This is supported by the DOS graphs shown in fig 6.11.
On the SXRD reported 2/3ML coverage surface, there is little total difference between the model shown from computational results and the SXRD-derived surface. The calculated total energy difference between these two surfaces is quite small, the computationally derived surface being only ~0.2eV below the SXRD-derived surface. However, when the SXRD model is subjected to simulated annealing, rather than just an energy minimisation, it reconstructs considerably, suggesting that the surface model is not stable. The computationally-derived model shown in fig 6.13b does reconstruct significantly on annealing.
The results show that the potassium atoms on the surface are of lower total energy when they are in differing rows and columns to one another, this is clear from figure 6.8, implying that the surface deposition of potassium builds up not in dimers, as proposed in ref , but in a more distributed way.
There is clear evidence that the 1ML surface requires more careful study, and that the surface needs to be probed in such a way as to be determine for certain whether the potassium layer or the silicon carbide underlayer is being examined by the STM results reported by Derycke et al..
The application of metallic overlayers on semiconductors are of paramount importance in the semiconductor industry. This study builds on the fundamentals of such studies. The skills developed in the study of the K-SiC(100) surface can be applied to the In-Si(111) system discussed in Chapter 8 as well as further studies of metals on the SiC(100) surface.
In the study of the 3C-(111)-SiC and the 6H-(0001)-SiC √3 × √3 surfaces, the band structure modifications over time show compelling evidence that the surface transition may be due to a soft phonon interaction. Previously this has been ascribed to a Mott-Hubbard transition. The evidence for soft phonons being the mechanism involved is not solely based on the band structure. The WHAM graphs in figures 7.13 and 7.14, although not quantitatively ideal, indicate two local minima on each of the surfaces, indicative of two metastable surface reconstructions and providing further evidence for the soft phonon hypothesis.
Experimentally, it is very difficult to find evidence to show the soft phonon theory as being the mechanism that drives this system. Although the computational evidence presented here supports published band structures, the metastable states identified require a more accurate calculation of the electronic structure for direct comparison with both Mott-Hubbard model calculated results and experimental results. During the course of the literature review on this surface, no optical data was found. If such studies were carried out both experimentally and theoretically the comparison may well yield a far greater understanding of the mechanisms involved here.
There are a number of surfaces that have been discussed in the literature as exhibiting a Mott-Hubbard transition. The Mott-Hubbard is an important and exciting phenomena in it’s own right, but it is not a catch-all solution, and much new and novel physics can be missed by categorising many different phenomena with the same brush. This study shows that in some cases the Mott-Hubbard transition is not the underlying mechanism and that what is actually occurring is elegant and an interesting problem to study.
The FIREBALL code has evolved considerably into FIREBALL-Lightning and Create-Lightning during the course of this work. The principle aim was to make it far easier to develop new algorithms and implement new functionality. This has been achieved successfully with the modularisation of the code. Rewriting the codebase was a substantial undertaking but worthwhile in that there are far fewer barriers to the adoption of FIREBALL within academic communities. Informal discussion with groups in other universities may result in the use of the Lightning code in massively parallel energy minimisation functions, for example using point energy to calculate a structure with no initial estimation of the cell structure or in using genetic algorithms for material design, using FIREBALL-Lightning for evaluation of each generation.Evaluation of the new codebase is ongoing at time of writing in addition to new applications of the code being concept tested.
The new applications that FIREBALL-Lightning can be applied to are endless. At the time of writing, work has begun to develop crowd-sourced heuristic techniques for material prediction that will incorporate the new, faster FIREBALL-Lightning as its primary computational tool. As well as this, there is work ongoing in computing an optical spectrum by the dipole summing technique and by the use of electrostatic potentials- both of which are possible in the new modular approach to the FIREBALL code.
Regarding the study of In-Si(111), this surface has enjoyed a lot of interest for the past decade or so, and its behaviour continues to be explored. The soft phonon and surface metallisation properties of this surface make it an ideal candidate for progression of the work described on silicon carbide in this thesis. There has been much work published on the optical properties of this system. There are theoretical and experimental RA spectra in the literature that can be used for direct comparison of a FIREBALL implemented calculation of this property. The completion of a FIREBALL implemented calculation would be a logical step forward from the work in this thesis as it would combines phonon interactions and surface metallisation with optical properties, all of which are discussed in this work.
Work is ongoing towards this end. The improvements of Create-Lightning and FIREBALL-Lightning packages are continually delivering speed increases over the older codebase by a factor of at least 2 and have spurned the effort considerably.
This Appendix merely consists of notes drawn up during the testing and implementation of a numerical quadrature method- adaptive Simpson’s algorithm into FIREBALL.
Adaptive Simpson’s method is based on estimating the error seen in calculating a definite integral using Simpson’s rule. If the error exceeds a user-specified tolerance, the algorithm calls for subdividing the interval of integration in two and applying adaptive Simpson’s method to each subinterval in a recursive manner.
This continues until:
Where:- a and b are the ends of an interval with midpoint c, S() is the Simpson’s rule estimate of that interval and ε is the relative error. The 15 ensures that estimates obtained are exact for polynomials of degree 5 or less. (This is a proof from ref [])
Essentially, if the error is more than a specific tolerance, the program halves the section, and then goes again, if this is more then the error, it splits the section again and repeats the process.
The following set of points is a pseudocode of the Adapative Simpson’s Algorithm, the flow-chart of which is depicted in figure 5.2.
a. Get upper and lower limits from calling program
b. Get error level that’s acceptable from calling program
c. Get function to be integrated from calling program
d. Evaluate integral limitations.
e. Evaluate first estimate
1. Simpsons rule applied as if there is only two sections
2. Check against Rule applied against an evaluation of just one section.
3. If difference is less than 15epsilon, done.
f. Evaluate second estimate
1. Left half function from last step split into two sections
2. checked against Rule applied against an evaluation of just one section.
3. If difference is less than 15epsilon, done.
4. Otherwise halve current section and repeat (6) until
g. Result is sum of all sections once barrier is reached of 15ε.
Table A.I.1 shows a range of input and output results for various functions used to test the Adaptive Simpson’s module. The results were compared to the result from the analytical solution of each integrand.
The deviation between analytical results and the program output can be seen in table A.I.1, the module results break down at polynomials of order 6 or so.
The derivation of exact exchange for FIREBALL is included in this appendix for completeness. The bulk of this work was carried out by Prof. J. Lewis (co-supervisor of this thesis) during a sabbatical in 2000.
The single-particle Hartree-Fock equation can be written as:
where hi is the single-particle Hamiltonian given by:
and V ex(r,r’) is the exchange-potential;
The exchange potential operating on orbital ψi(ri) can be taken as (for our purposes)
Assume that the wavefunctions ψ are expanded as a linear combination of atomic orbitals, such that ψj(r) = ∑ cjμ ψμ(r - r’). Therefore in the LCAO formalism the exchange potential operating on an atomic orbital will be
We are interested in the expectation value between two atomic orbitals < ψμ (r - r1)|Vex ν(r)| ψν (r - r2)> , which is:
Using relation (3), at the end of this appendix, for the expansion of and noting that each orbital can be separated into a radial piece and an angular piece. For 1C r1 can be set to zero.
The angular integrals can be evaluated based on relations (1) and (2), at the end of this appendix.
Multiplication of the two angular integrals yields:
It is confirmed from the Clebsch-Gordan coefficients that only the μ=υ terms survive. This is what we would expect. Therefore, the angular terms reduce to:
The net result is then:
The two-center integrals for the ontop (L) case are evaluated as follows (rewritten from previous equations):
Use relation (3) for the expansion of and note that each orbital can be separated into a radial piece and an angular piece. Also, change variable such that r → r’ and r’ → r. Note: r1 = 0 and r2 = dz can be used in this case.
The angular integral for can be evaluated based on relations (1) and (2)
The radial integral for r’ is
Now consider the remaining integral and note that
The ψ integration reduces to
The remaining integral can be rewritten as an integral over ρ and z such as has been done for most two-center integrals
For example, in a p-state,
The two-center integrals for the ontop (R) case are evaluated as follows:
Using relation (3) for the expansion of and note that each orbital can be separated into a radial piece and angular piece. Note that r1 = 0 and r2 = dz is used for this case
The angular integral for can be evaluated based on relations (1) and (2)
The radial integral for r’ is
Consider the remaining integral and note that
The ψ integration reduces to
The remaining integral can be rewritten as an integral over ρ and z as for the ontop (L) case,
The atom case unfortunately cannot be analytically simplified as well as the ontop cases. For this case, we have:
Using relation (3) for the expansion of , this becomes
The only integral which can be used analytically is the integration (cylindrical coordinates) we have
The remaining integrals are done in cylindrical coordinates, yielding the two resulting integrals to be done numerically.
Now consider the case where α = β only, and r3 = r4. This is a case where essentially the potential V Xv (r, r’) is taken for wavefunctions centered at the same site. These would be the “diagonal” terms in the potential. For this case the matrix element reduces to
This should be the largest three-center interaction. Currently, we are only interested in one- and two- center interactions, so reduce the matrix element equation further for these special cases. Later maybe consider the case where α ≠ β (but r3 = r4)
(L) Ontop Case r3 = r4
(R) Ontop Case r3 = r2
Atom Case r1 = r2
Relating angle to
(3) Expansion at
(see C.Cohen Tannoudjij1 p1046)
A specific equation for the Clebsch-Gordon coefficients is given in ref. , page 178.
Note that the l values are bound by |l1 - l2| ≤ l ≤ l1 + l2 . Possible values of l1, l2 and l have been previously worked out for kinetic matrix elements. Note that z spans all integers (positive and negative), but factorial involving z must be non-negative in order for the corresponding term in the sum to be included.
Here is the table of possible l values which survive; based on the fact that l1 = 0,1,2,3 and l2 = 0,1,2,3 which corresponds to s, p, d, and f orbitals.
With the condition that any factorial involving z must be non negative, determine the allowable values of z given the allowable values of l1m1, l2m2, and lm. The following conditions must be upheld:
The recursion relations for the Clebsch-Gordon coefficients are as follows:
To generate the associated Legendre Polynomials
 C. Cohen-Tannoudji, B. Diu and F. Laloe, “Quantum Mechanics”, Wiley-Interscience, 2006, ISBN: 9780471569527.
 D. Bohm, “Quantum Theory”, New York: Prentice Hall, reprint, 1989, ISBN: 048659690
In order to calculate the binding energy of the potassium to the surface in chapter 6, it was necessary to calculate the charge migration per potassium atom. In figure A.III.1, the average charge lost per potassium atom with standard deviation errorbars is shown for a number of studies for each coverage simulated.
As mentioned in Chapter 6, the z-axis position of the potassium above the silicon carbide surface was also calculated, figure A.III.2 shows the calculated average distance in the z direction of the potassium atoms for a number of simulations per coverage. The errorbars show the standard deviations from the average.
In figure 6.9, a scatter plot showing the annealed surfaces cohesive energy versus coverage. This neglects to show the energies calculated for the non-annealed (quenched only) surfaces, figure A.III.3 shows both the non annealed and annealed results for coverage versus cohesive energy.
Figure 6.10 shows the surface configurations that were of lowest energy and could be compared to experimental results. If figure A.III.4, the chemical potential plot for potassium is shown for many more annealed surface reconstructions. In the legend for figure A.III.4, each colour is defined by the following labeling system: the number is the number of atoms on the 12 pedestal-site silicon carbide surface, if there is no letter after the number, the potassium atoms are randomly-placed, a letter “a” indicated the potassium was placed in the direction of chains reported by Derycke et al., a “b” indicates the atoms are placed so as to form chains perpendicular on the surface to “a”, and “c” indicates that atoms were placed so as to form the asymmetric dimer rows seen in figure 6.13 for the 1ML coverage (6.13(a)).
Many more DOS plots were generated than are reported on in this thesis, those of note that did not ideally fit the flow of Section 6.2.3 are shown in figure A.III.5, including surface DOS of coverages of greater than 1ML. In figure A.III.5 the semiconducting DOS of the silicon carbide in the bulk is also shown.
In figures A.III.6 - A.III.8 the lowest energy surface reconstructions found for coverages greater than 1ML is shown, a side view as well as top view is shown to express the ejection of excess potassium atoms on the surface.
Figure 7.7 show a summed band structure over a single pass of a soft phonon for the 2×(√3× √3) surface, which has four surface adatoms. Figures A.IV.1-A.IV.2 are the band structures calculated in the same way for the √3× √3 surface.
In another set of experiments, the supercell consisted of a 2√3×√3 and a √3×2√3 surface. Figures A.IV.3 and A.IV.4 show each of the summed band structures for these studies.
As with figure 7.7, these band structures show a band within the gap that moves with time, resulting in a smearing of that band, which may explain the splitting of the band seen experimentally.
During this project, and to help with the rewriting of the Create code into Create-Lightning, it was decided to make a flow chart of the 2006 code by meticulously going through the entire package. This is the flowchart of Create-2006, with branches for some, but not all, of the “if” statements within the code. By comparison a flow chart of the Create-Lightning code can be seen in figure 5.1.
 F. Pedreschi, J. D. O’Mahony, P. Weightman, J. R. Power. Evidence of electron confinement in the single-domain (4 x 1)-in superstructure on vicinal si(111). Applied Physics Letters, 73(15):2152–2154, 1998.
 Vincent Derycke, Patrick G. Soukiassian, Fabrice Amy, Yves J. Chabal, Marie D. D’angelo, Hanna B. Enriquez, Mathieu G. Silly. Nanochemistry at the atomic scale revealed in hydrogen-induced semiconductor surface metallization. Nat Mater, 2(4):253–258, 04 2003.
 A. P. Favaro, K. Capelle, Joao Batista Ferreira. Construction of model dielectric functions for two- and three-dimensional electron liquids from density functionals. Phys. Rev. B, 73(4):045133, Jan 2006.
 S. F. Boys. Electronic wave functions. i. a general method of calculation for the stationary states of any molecular system. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 200(1063):542–554, 1950.
 M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias, J. D. Joannopoulos. Iterative minimization techniques for ab initio total-energy calculations: molecular dynamics and conjugate gradients. Rev. Mod. Phys., 64(4):1045–1097, Oct 1992.
 M D Segall, Philip J D Lindan, M J Probert, C J Pickard, P J Hasnip, S J Clark, M C Payne. First-principles simulation: ideas, illustrations and the castep code. Journal of Physics: Condensed Matter, 14(11):2717, 2002.
 X Gonze, G Rignanese, M Verstraete, J Betiken, Y Pouillon, R Caracas, F Jollet, M Torrent, G Zerah, M Mikami, et al. A brief introduction to the abinit software package. Zeitschrift fur Kristallographie, 220(5-6-2005):558–562, 2005.
 Jose M Soler, Emilio Artacho, Julian D Gale, Alberto Garcia, Javier Junquera, Pablo Ordejen, Daniel Sanchez-Portal. The siesta method for ab initio order- n materials simulation. Journal of Physics: Condensed Matter, 14(11):2745, 2002.
 Otto F. Sankey, David J. Niklewski. Ab initio multicenter tight-binding model for molecular-dynamics simulations and other applications in covalent systems. Phys. Rev. B, 40(6):3979–3995, Aug 1989.
 James P. Lewis, Kurt R. Glaesemann, Gregory A. Voth, Jurgen Fritsch, Alexander A. Demkov, Jose Ortega, Otto F. Sankey. Further developments in the local-orbital density-functional-theory tight-binding method. Phys. Rev. B, 64(19):195103, Oct 2001.
 Pavel Jelinek, Hao Wang, James P. Lewis, Otto F. Sankey, José Ortega. Multicenter approach to the exchange-correlation interactions in ab initio tight-binding methods. Phys. Rev. B, 71(23):235101, Jun 2005.
 Otto F. Sankey, Alexander A. Demkov, Wolfgang Windl, JÃŒrgen H. Fritsch, James P. Lewis, Miguel Fuentes-Cabrera. The application of approximate density functionals to complex systems. International Journal of Quantum Chemistry, 69(3):327–340, 1998.
 X. Gonze, M. Scheffler. Exchange and correlation kernels at the resonance frequency: Implications for excitation energies in density-functional theory. Phys. Rev. Lett., 82(22):4416–4419, May 1999.
 Patrick G Soukiassian, Hanna B Enriquez. Atomic scale control and understanding of cubic silicon carbide surface reconstructions, nanostructures and nanochemistry. Journal of Physics: Condensed Matter, 16(17):S1611, 2004.
 F. Bechstedt, P. Kackell, A. Zywietz, K. Karch, B. Adolph, K. Tenelsen, J. Furthmüller. Polytypism and properties of silicon carbide. Physica Status Solidi B Basic Research, 202:35–62, Lipiec 1997.
 P. Soukiassian, F. Semond, L. Douillard, A. Mayne, G. Dujardin, L. Pizzagalli, C. Joachim. Direct observation of a beta-sic(100)- c(4x2) surface reconstruction. Phys. Rev. Lett., 78(5):907–910, Feb 1997.
 V. Derycke, P. Fonteneau, Y. K. Hwu, P. Soukiassian. From k atom pairs to k atomic chains: A semiconducting 2 x 3 to metallic 2 x 1 transition on the beta-sic(100) c(4 x 2) surface. Applied Physics Letters, 88(2):022105, 2006.
 Alessandra Catellani, Giulia Galli, Francois Gygi, Fabio Pellacini. Influence of stress and defects on the silicon-terminated sic(001) surface structure. Phys. Rev. B, 57(19):12255–12261, May 1998.
 V. Yu. Aristov, P. Soukiassian, A. Catellani, R. Di Felice, G. Galli. Experimental and theoretical electronic structure determination of the beta-sic(001)c(4x2) surface reconstruction. Phys. Rev. B, 69(24):245326, Jun 2004.
 Fernando B Mota, Von B Nascimento, Caio M C de Castilho. Ab initio electronic and structural properties of clean and hydrogen saturated sic(100)(3 x 2) surfaces. Journal of Physics: Condensed Matter, 17(30):4739, 2005.
 R. Di Felice, C. M. Bertoni, C. A. Pignedoli, A. Catellani. Hydrogen-induced surface metallization of β - sic(100) - (3 × 2) revisited by density functional theory calculations. Phys. Rev. Lett., 94(11):116103, Mar 2005.
 A. Franciosi, P. Philip, S. Chang, A. Wall, A. Raisanen, N. Troullier, P. Soukiassian. Electronic promoters and semiconductor oxidation: Alkali metals on si(111) surfaces. Phys. Rev. B, 35(2):910–913, Jan 1987.
 A. Tejeda, E. Wimmer, P. Soukiassian, D. Dunham, E. Rotenberg, J. D. Denlinger, E. G. Michel. Atomic structure determination of the 3c-sic(001) c(4x2) surface reconstruction: Experiment and theory. Phys. Rev. B, 75(19):195315, May 2007.
 Peter Deak, Balint Aradi, Jan M. Knaup, Thomas Frauenheim. Hydrogen adsorption and etching on the si-rich 3 c -sic(001) 3 × 2 surface: First-principles molecular dynamics calculations. Phys. Rev. B, 79(8):085314, Feb 2009.
 Xiangyang Peng, Peter Kruger, Johannes Pollmann. Metallization of the 3c-sic(001)-(3x2) surface induced by hydrogen adsorption: A first-principles investigation. Phys. Rev. B, 72(24):245320, Dec 2005.
 L. Li, Y. Hasegawa, T. Sakurai. Si- and c-rich structure of the 6h-sic(0001) surface. Journal of Vacuum Science and Technology B: Microelectronics and Nanometer Structures, 15(4):1307–1309, 1997. cited By (since 1996) 8.
 J. Furthmuller, P. Kackell, F. Bechstedt, G. Kresse. Extreme softening of vanderbilt pseudopotentials: General rules and case studies of first-row and d-electron elements. Phys. Rev. B, 61(7):4576–4587, Feb 2000.
 P. Kackell, J. FurthmÂžller, F. Bechstedt. Polytypism and surface structure of sic. Diamond and Related Materials, 6(10):1346 – 1348, 1997. Proceeding of the 1st European Conference on Silicon Carbide and Related Materials (ECSCRM 1996).
 L. Li, Y. Hasegawa, T. Sakurai, I. S. T. Tsong. Field-ion scanning tunneling microscopy study of the atomic structure of 6h–sic(0001) surfaces cleaned by in situ si molecular beam etching. Journal of Applied Physics, 80(4):2524–2526, 1996.
 L. S. O. Johansson, L. Duda, M. Laurenzis, M. Krieftewirth, B. Reihl. Electronic structure of the 6h-sic(0001)-3x3 surface studied with angle-resolved inverse and direct photoemission. Surface Science, 445(1):109 – 114, 2000.
 U. Starke, J. Schardt, J. Bernhardt, M. Franke, K. Reuter, H. Wedler, K. Heinz, J. Furthmuller, P. Kackell, F. Bechstedt. Novel reconstruction mechanism for dangling-bond minimization: Combined method surface structure determination of sic(111)- (3x3). Phys. Rev. Lett., 80(4):758–761, Jan 1998.
 J. Hubbard. Electron correlations in narrow energy bands. iii. an improved solution. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 281(1386):401–419, 1964.
 F. Gebhard. The Mott Metal-Insulator Transition: Models and Methods. Springer, Germany, 1997. ISBN: 9783540614814@ArticlePhysRevLett.84.135, title = U Parameter of the Mott-Hubbard Insulator 6H -SiC(0001)-(×)R30∘: An Ab Initio Calculation, author = Rohlfing, Michael and Pollmann, Johannes , journal = Phys. Rev. Lett., volume = 84, number = 1, pages = 135–138, numpages = 3, year = 2000, month = Jan, doi = 10.1103/PhysRevLett.84.135, publisher = American Physical Society .
 F. Flores, J. Ortega, R. PÃⒸrez. Many-body effects and the metal-insulator transition at semiconductor surfaces and interfaces. Surface Review and Letters, 6(3-4):411–433, 1999. cited By (since 1996) 22.
 D. FarÃas, W. Kaminski, J. Lobo, J. Ortega, E. Hulpke, R. PÃⒸrez, F. Flores, E. G. Michel. Phonon softening, chaotic motion, and order-disorder transition in sn/ge(111). Phys. Rev. Lett., 91(1):016103, Jul 2003.
 J. Avila, A. Mascaraque, E. G. Michel, M. C. Asensio, G. LeLay, J. Ortega, R. PÃⒸrez, F. Flores. Dynamical fluctuations as the origin of a surface phase transition in sn/ge(111). Phys. Rev. Lett., 82(2):442–445, Jan 1999.
 C GonzÃ¡lez, J Ortega, F Flores. Metal-insulator transition in one-dimensional in-chains on si(111): combination of a soft shear distortion and a double-band peierls instability. New Journal of Physics, 7(1):100, 2005.
 Hong Wang, Jessica Stalnaker, Hilaire Chevreau, James P. Lewis. Potential of mean force calculations using ab initio tight-binding molecular dynamics: Application to n-no2 bond dissociation in dmna and hmx. Chemical Physics Letters, 457(1-3):26 – 30, 2008.
 H. W. Yeom, S. Takeda, E. Rotenberg, I. Matsuda, K. Horikoshi, J. Schaefer, C. M. Lee, S. D. Kevan, T. Ohta, T. Nagao, S. Hasegawa. Instability and charge density wave of metallic quantum chains on a silicon surface. Phys. Rev. Lett., 82(24):4898–4901, Jun 1999.
 A. A. Saranin, A.V. Zotov, K. V. Ignatovich, V. G. Lifshits, T. Numata, O. Kubo, H. Tani, M. Katayama, K. Oura. Structural model for the si(111)-4x1-in reconstruction. Phys. Rev. B, 56(3):1017–1020, Jul 1997.
 A. A. Saranin, A. V. Zotov, V. G. Lifshits, J. T. Ryu, O. Kubo, H. Tani, T. Harada, M. Katayama, K. Oura. Analysis of surface structures through determination of their composition using stm: Si(100)4x3-in and si(111)4x1-in reconstructions. Phys. Rev. B, 60(20):14372–14381, Nov 1999.
 C. Kumpf, O. Bunk, J. H. Zeysing, Y. Su, M. Nielsen, R. L. Johnson, R. Feidenhans’l, K. Bechgaard. Low-temperature structure of indium quantum chains on silicon. Phys. Rev. Lett., 85(23):4916–4919, Dec 2000.
 O. Bunk, G. Falkenberg, L. Seehofer, J. H. Zeysing, R. L. Johnson, M. Nielsen, R. Feidenhans’l, E. Landemark. Structure determination of the indium induced si(001)-(4 x 3) reconstruction by surface x-ray diffraction and scanning tunneling microscopy. Applied Surface Science, 123-124:104 – 110, 1998. Proceedings of the Sixth International Conference on the Formation of Semiconductor Interfaces.
 O. Bunk, G. Falkenberg, J. H. Zeysing, L. Lottermoser, R. L. Johnson, M. Nielsen, F. Berg-Rasmussen, J. Baker, R. Feidenhans’l. Structure determination of the indium-induced si(111)-(4x1) reconstruction by surface x-ray diffraction. Phys. Rev. B, 59(19):12228–12231, May 1999.
 A. A. Saranin, E. A. Khramtsova, K. V. Ignatovich, V. G. Lifshits, T. Numata, O. Kubo, M. Katayama, I. Katayama, K. Oura. Indium-induced si(111)4x1 silicon substrate atom reconstruction. Phys. Rev. B, 55(8):5353–5359, Feb 1997.
 Helmut Ofner, Svetlozar L. Surnev, Yoram Shapira, Falko P. Netzer. Spectroscopy of interface states of indium-si(111)(4 x 1) and (1 x 1)r30 surfaces. Surf. Sci., 307-309(Part 1):315 – 320, 1994.
 M. K. Kelly, G. Margaritondo, J. Anderson, D. J. Frankel, G. J. Lapeyre. Reconstruction of aluminum and indium overlayers on si(111): A systematic study with high-resolution electron energy loss spectroscopy and low-energy electron diffraction. Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films, 4(3):1396–1399, 1986.
 T. Abukawa, M. Sasaki, F. Hisamatsu, T. Goto, T. Kinoshita, A. Kakizaki, S. Kono. Surface electronic structure of a single-domain si(111)4 x 1-in surface: a synchrotron radiation photoemission study. Surface Science, 325(1-2):33 – 44, 1995.
 Kazuyuki Sakamoto, Hidenori Ashima, Han Woong Yeom, Wakio Uchida. Angle-resolved high-resolution electron-energy-loss study of in-adsorbed si(111)-(4x1) and -(8x2) surfaces. Phys. Rev. B, 62(15):9923–9926, Oct 2000.
 S. Chandola, K. Hinrichs, M. Gensch, N. Esser, S. Wippermann, W. G. Schmidt, F. Bechstedt, K. Fleischer, J. F. McGilp. Structure of si(111)-in nanowires determined from the midinfrared optical response. Phys. Rev. Lett., 102(22):226805, Jun 2009.
 C. Gonzalez, Jiandong Guo, J. Ortega, F. Flores, H. H. Weitering. Mechanism of the band gap opening across the order-disorder transition of si(111)(4x1)-in. Phys. Rev. Lett., 102(11):115501, Mar 2009.
 S. Wippermann, W.G. Schmidt, A. Calzolari, M. Buongiorno Nardelli, A.A. Stekolnikov, K. Seino, F. Bechstedt. Quantum conductance of in nanowires on si(111) from first principles calculations. Surface Science, 601(18):4045 – 4047, 2007. ECOSS-24, Proceedings of the 24th European Conference on Surface Science.
 D. Ronnow, L. F. Lastras-Martinez, M. Cardona, P. V. Santos. Determination of the piezo-optical properties of semiconductors above the fundamental gap by means of reflectance difference spectroscopy. J. Opt. Soc. Am. A, 16(3):568–573, Mar 1999.
 S. Selci, A. Cricenti, F. Ciccacci, A.C. Felici, C. Goletti, Zhu Yong, G. Chiarotti. Dielectric functions of si(111)2x1, ge(111)2x1, gaas(110) and gap(110) surfaces obtained by polarized surface differential reflectivity. Surface Science, 189-190:1023 – 1027, 1987. Proceedings of the Ninth European Conference on Surface Science.
 R. J. Cole, B. G. Frederick, P. Weightman. Substrate dependence of adlayer optical response in reflectance anisotropy spectroscopy. Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films, 16(5):3088–3095, 1998.
 K. Hingerl, D. E. Aspnes, I. Kamiya, L. T. Florez. Relationship among reflectance-difference spectroscopy, surface photoabsorption, and spectroellipsometry. Applied Physics Letters, 63(7):885–887, 1993.
 D. E. Aspnes. Above-bandgap optical anisotropies in cubic semiconductors: A visible–near ultraviolet probe of surfaces. Journal of Vacuum Science & Technology B: Microelectronics and Nanometer Structures, 3(5):1498–1506, 1985.
 M. Wassermeier, I. Kamiya, D. E. Aspnes, L. T. Florez, J. P. Harbison, P. M. Petroff. Optical spectroscopy of (001) gaas and alas under molecular-beam epitaxy growth conditions. Journal of Vacuum Science & Technology B: Microelectronics and Nanometer Structures, 9(4):2263–2267, 1991.
 Patrizia Monachesi, Maurizia Palummo, Rodolfo Del Sole, Rajeev Ahuja, Olle Eriksson. Reflectance anisotropy spectra of cu and ag (110) surfaces from ab initio theory. Phys. Rev. B, 64(11):115421, Aug 2001.
 Wenchang Lu, W. G. Schmidt, E. L. Briggs, J. Bernholc. Optical anisotropy of the sic(001)- (3x2) surface: Evidence for the two-adlayer asymmetric-dimer model. Phys. Rev. Lett., 85(20):4381–4384, Nov 2000.
 Conor Hogan, Rodolfo Del Sole, Giovanni Onida. Optical properties of real surfaces from microscopic calculations of the dielectric function of finite atomic slabs. Phys. Rev. B, 68(3):035405, Jul 2003.
 W. G. Schmidt, N. Esser, A. M. Frisch, P. Vogt, J. Bernholc, F. Bechstedt, M. Zorn, Th. Hannappel, S. Visbeck, F. Willig, W. Richter. Understanding reflectance anisotropy: Surface-state signatures and bulk-related features in the optical spectrum of inp(001)(2 × 4). Phys. Rev. B, 61(24):R16335–R16338, Jun 2000.
 M. Makowska-Janusik, J. Sanetra, H. Palmer, D. Bogdal, E. Gondek, I. V. Kityk. Absorption spectra of poly-n-vinylcarbazole derivatives by experiment and simulation. European Polymer Journal, 40(4):799 – 804, 2004.