Creation of Ordered Layers on Semiconductor Surfaces: An ab Initio Molecular Dynamics Study of the SiC(001)-3×2 and SiC(100)-c(2×2) Surfaces

The chemistry of hybrid structures composed of organic molecules and semiconductor surfaces is opening up exciting new areas of development in molecular electronics, nanoscale sensing devices, and surface lithography (Filler & Bent, 2002; Kachian et al., 2010; Kruse & Wolkow, 2002). Covalent attachment of organic molecules to a semiconducting surface can yield active devices, such as molecular switches (Filler & Bent, 2003; Flatt et al., 2005; Guisinger, Basu, Greene, Baluch & Hersam, 2004; Guisinger, Greene, Basu, Baluch & Hersam, 2004; He et al., 2006; Rakshit et al., 2004) and sensors (Cattaruzza et al., 2006) or passivating insulating layers. Moreover, it is assumed that the reactions can be controlled by “engineering” specific modifications to organic molecules, suggesting possible new lithographic techniques. One of the goals of controlling the surface chemistry is the creation of ordered nanostructures on semiconducting surfaces. Indeed, there has been some success in obtaining locally ordered structures on the hydrogen terminated Si(100) surface (Basu et al., 2006; Hossain et al., 2005b; Kirczenow et al., 2005; Lopinski et al., 2000; Pitters et al., 2006). These methods require a dangling Si bond without a hydrogen to initialize the self-replicating reaction. Another popular approach eliminates the initialization step by exploiting the reactivity between surface dimers on certain reconstructed surfaces with the π bonds in many organic molecules. The challenge with this approach lies in designing the surface and/or the molecule so as to eliminate all but one desired reaction channel. Charge asymmetries, such as occur on the Si(100)-2×1 surface, lead to a violation of the usual Woodward-Hoffman selection rules, which govern many purely organic reactions, allowing a variety of possible [4+2] and [2+2] surface adducts. Silicon-carbide (SiC) is often the material of choice for electronic and sensor applications under extreme conditions (Capano & Trew, 1997; Mélinon et al., 2007; Starke, 2004) or subject to biocompatibility constraints (Stutzmann et al., 2006). SiC has a variety of reconstructions that could possibly serve as candidates for creating ordered organic/semiconductor interfaces. However, the choice of the reconstruction is crucial. Multiple reactive sites, such as occur on some of the SiC surfaces, will lead to a broad distribution of adducts. The exploration of hybrid organic-semiconductor materials and the reactions associated with them is an area in which theoretical and computational tools can play an important role. Indeed, modern theoretical methods combined with high-performance computing, 10


Introduction
The chemistry of hybrid structures composed of organic molecules and semiconductor surfaces is opening up exciting new areas of development in molecular electronics, nanoscale sensing devices, and surface lithography (Filler & Bent, 2002;Kachian et al., 2010;Kruse & Wolkow, 2002).
Covalent attachment of organic molecules to a semiconducting surface can yield active devices, such as molecular switches (Filler & Bent, 2003;Flatt et al., 2005;Guisinger, Basu, Greene, Baluch & Hersam, 2004;Guisinger, Greene, Basu, Baluch & Hersam, 2004;He et al., 2006;Rakshit et al., 2004) and sensors (Cattaruzza et al., 2006) or passivating insulating layers.Moreover, it is assumed that the reactions can be controlled by "engineering" specific modifications to organic molecules, suggesting possible new lithographic techniques.One of the goals of controlling the surface chemistry is the creation of ordered nanostructures on semiconducting surfaces.Indeed, there has been some success in obtaining locally ordered structures on the hydrogen terminated Si(100) surface (Basu et al., 2006;Hossain et al., 2005b;Kirczenow et al., 2005;Lopinski et al., 2000;Pitters et al., 2006).These methods require a dangling Si bond without a hydrogen to initialize the self-replicating reaction.Another popular approach eliminates the initialization step by exploiting the reactivity between surface dimers on certain reconstructed surfaces with the π bonds in many organic molecules.The challenge with this approach lies in designing the surface and/or the molecule so as to eliminate all but one desired reaction channel.Charge asymmetries, such as occur on the Si(100)-2×1 surface, lead to a violation of the usual Woodward-Hoffman selection rules, which govern many purely organic reactions, allowing a variety of possible [4+2] and [2+2] surface adducts.Silicon-carbide (SiC) is often the material of choice for electronic and sensor applications under extreme conditions (Capano & Trew, 1997;Mélinon et al., 2007;Starke, 2004) or subject to biocompatibility constraints (Stutzmann et al., 2006).SiC has a variety of reconstructions that could possibly serve as candidates for creating ordered organic/semiconductor interfaces.However, the choice of the reconstruction is crucial.Multiple reactive sites, such as occur on some of the SiC surfaces, will lead to a broad distribution of adducts.The exploration of hybrid organic-semiconductor materials and the reactions associated with them is an area in which theoretical and computational tools can play an important role.Indeed, modern theoretical methods combined with high-performance computing, have advanced to a level such that the thermodynamics and reaction mechanisms can be routinely studied.These studies can aid in the interpretation of experimental results and can leverage theoretical mechanisms to predict the outcomes of new experiments.This chapter will focus on a description of one set of such techniques, namely, those based on density functional theory and first-principles or ab initio molecular dynamics (Car & Parrinello, 1985;Marx & Hutter, 2000;2009;Tuckerman, 2002).As these methods employ an explicit representation of the electronic structure, electron localization techniques can be used to follow local electronic rearrangements during a reaction and, therefore, generate a clear picture of the reaction mechanism.In addition, statistical mechanical tools can be employed to obtain thermodynamic properties of the reaction products, including relative free energies and populations of the various products.Following a detailed description of the computational approaches, we will present two applications of conjugated dienes reacting with different choices of SiC surfaces (Hayes & Tuckerman, 2008).We will investigate how the surface structure influences the thermodynamics of the reaction products and how these thermodynamic properties can be used to guide the choice of the surface in order to control the product distribution and associated free energies.

Computational methods
Because the problem of covalently attaching an organic molecule to a semiconductor surface requires the formation of chemical bonds, a theoretical treatment of this problem must be able to describe this bond formation process, which generally requires an ab initio approach in which the electronic structure is accounted for explicitly.Assuming the validity of the Born-Oppenheimer approximation, the goal of any ab initio approach is to approximate the ground-state solution of the electronic Schrödinger equation Ĥelec (R)|Ψ 0 (R) = E 0 (R)|Ψ 0 (R) ,w h e r e Ĥelec is the electronic Hamiltonian, and R denotes a classical configuration of the chemical nuclei in the system.Ultimately, in order to predict reaction mechanisms and thermodynamics, one needs to use the approximate solution of the Schrödinger equation to propagate the nuclei dynamically using Newton's laws of motion.This is the essence of the method known as ab initio molecular dynamics (AIMD) (Car & Parrinello, 1985;Marx & Hutter, 2000;2009;Tuckerman, 2002).Unless otherwise stated, all of the calculations to be presented in this chapter were carried out using the implementation of plane-wave based AIMD implemented in the PINY_MD package (Tuckerman et al., 2000).In this section, we will briefly review density functional theory (DFT) as the electronic structure method of choice for the studies to be described in this chapter.DFT represents an optimal compromise between accuracy and computational efficiency.This is an important consideration, as a typical AIMD calculation requires that the electronic structure problem be solved tens to hundreds of thousands of time in order to generate one or more trajectories of sufficient length to extract dynamic and thermodynamic properties.We will then describe the AIMD approach, including several technical considerations such as basis sets, boundary conditions, and electron localization schemes.

Density functional theory
As noted above, we seek approximate solutions to the electronic Schrödinger equation.To this end, we begin by considering a system of N nuclei at positions R 1 , ..., R N ≡ R and M electrons with coordinate labels r 1 , ..., r M and spin states s 1 , ..., s M .The fixed nuclear positions allow us to define the electronic Hamiltonian (in atomic units) as where Z I is the charge on the Ith nucleus.The time-independent electronic Schrödinger equation or electronic eigenvalue problem which assumes the validity of the Born-Oppenheimer approximation, could, in principle, yield all of the electronic energy levels and eigenfunctions at the given nuclear configuration R.H e r e x i = r i , s i is a combination of coordinate and spin variables.Unfortunately, for large condensed-phase problems of the type to be considered here, an exact solution of the electronic eigenvalue problem is computationally intractable.The Kohn-Sham (KS) (Kohn & Sham, 1965) formulation of density functional theory (DFT) (Hohenberg & Kohn, 1964) replaces the fully interacting electronic system described by Eq. ( 2) by an equivalent non-interacting system that is required to yield the same ground-state energy and wave function as the original interacting system.As the name implies, the central quantity in DFT is the ground-state density n 0 (r) generated from the ground-state wave function Ψ 0 via where, for notational convenience, the dependence on R is left off.The central theorem of DFT is the Hohenberg-Kohn theorem, which states that there exists an exact energy E[n] that is a functional of electronic densities n(r) such that when E[n] is minimized with respect to n(r) subject to the constraint that dr n(r)=M (each n(r) must yield the correct number of electrons), the true ground-state density n 0 (r) is obtained.The true ground-state energy is then given by E 0 = E[n 0 ].The KS noninteracting system is constructed in terms of a set of mutually orthogonal single-particle orbitals ψ i (r) in terms of which the density n(r) is given by where f i are the occupation numbers of a set of N s such orbitals, where ∑ i f i = M.I n closed-shell systems, the orbitals are all doubly occupied so that N s = M/2, and In open-shell systems, we treat all of the electrons in double and singly occupied orbitals explicitly and take N s = M.When virtual or unoccupied orbitals are needed, we can take N s > M/2 or N s > M for closed and open-shell systems, respectively, and take f i = 0forthe virtual orbitals.
In KS theory, the energy functional is taken to be The first term in the functional represents the noninteracting quantum kinetic energy of the electrons, the second term is the direct Coulomb interaction between two charge distributions, the third therm is the exchange-correlation energy, whose exact form is unknown, and the fourth represents the "external" Coulomb potential on the electrons due to the fixed nuclei, V ext (r, R)=− ∑ I Z I /|r − R I |.Minimization of Eq. ( 5) with respect to the orbitals subject to the orthogonality constraint leads to a set of coupled self-consistent field equations of the form where the KS potential V KS (r) is given by and λ ij is a set of Lagrange multipliers used to enforce the orthogonality constraint ψ i |ψ j = δ ij .If we introduce a unitary transformation U that diagonalizes the matrix λ ij into Eq.( 6), then we obtain the Kohn-Sham equations in the form where φ i (r)=∑ j U ij ψ j (r) are the KS orbitals and ε i are the KS energy levels, i.e., the eigenvalues of the matrix λ ij .If the exact exchange-correlation functional were known, the KS theory would be exact.However, because E xc [n] is unknown, approximations must be introduced for this term in practice.The accuracy of DFT results depends critically on the quality of the approximation.One of the most widely used forms for E xc [n] is known as the generalized-gradient approximation (GGA), where in E xc [n] is approximated as a local functional of the form where the form of the function f GGA determines the specific GGA approximation.

Ab initio molecular dynamics
Solution of the KS equations yields the electronic structure at a set of fixed nuclear positions R 1 , ..., R N ≡ R. Thus, in order to follow the progress of a chemical reaction, we need an approach that allows us to propagate the nuclei in time.If we assume the nuclei can be treated as classical point particles, then we seek the nuclear positions R 1 (t), ..., R N (t) as functions of time, which are given by Newton's second law where M I and F I are the mass and total force on the Ith nucleus.If the exact ground-state wave function Ψ 0 (R) were known, then the forces would be given by the Hellman-Feynman theorem Within the framework of KS DFT, the force expression becomes The equations of motion, Eq. ( 10), are integrated numerically for a set of discrete times t = 0, Δt,2Δt, ..., N Δt subject to a set of initial coordinates R 1 (0), ..., R N (0) and velocities Ṙ1 (0), ..., ṘN (0) using a solver such as the velocity Verlet algorithm: where F I (0) and F I (Δt) are the forces at t = 0a n dt = Δt, respectively.Iteration of Eq. ( 14) yields a full trajectory of N steps.Eqs. ( 13) and ( 14) suggest an algorithm for generating the finite-temperature dynamics of a system using forces generated from electronic structure calculations performed "on the fly" as the simulation proceeds: Starting with the initial nuclear configuration, one minimizes the KS energy functional to obtain the ground-state density, and Eq. ( 13) is used to obtain the initial forces.These forces are then used to propagate the nuclear positions to the next time step using the first of Eqs. ( 14).At this new nuclear configuration, the KS functional is minimized again to obtain the new ground-state density and forces using Eq. ( 13), and these forces are used to propagate the velocities to time t = Δt.
These forces can also be used again to propagate the positions to time t = 2Δt.Th epr oc edu r e is iterated until a full trajectory is generated.This approach is known as "Born-Oppenheimer" dynamics because it employs, at each step, an electronic configuration that is fully quenched to the ground-state Born-Oppenheimer surface.
An alternative to Born-Oppenheimer dynamics is the Car-Parrinello (CP) method (Car & Parrinello, 1985;Marx & Hutter, 2000;Tuckerman, 2002).In this approach, an initially minimized electronic configuration is subsequently "propagated" from one nuclear configuration to the next using a fictitious Newtonian dynamics for the orbitals.In this "dynamics", the orbitals are given a small amount of thermal kinetic energy and are made "light" compared to the nuclei.Under these conditions, the orbitals actually generate a potential of mean force surface that is very close to the true Born-Oppenheimer surface.The equations of motion of the CP method are where μ is a mass-like parameter for the orbitals (which actually has units of energy× time 2 ), and λ ij is the Lagrange multiplier matrix that enforces the orthogonality of the orbitals as a holonomic constraint on the fictitious orbital dynamics.Choosing μ small ensures that the orbital dynamics is adiabatically decoupled from the true nuclear dynamics, thereby allowing the orbitals to generate the aforementioned potential of mean force surface.For a detailed analysis of the CP dynamics, see Marx et al. (1999); Tuckerman (2002).As an illustration of the CP dynamics, Fig. 1 of Tuckerman & Parrinello (1994) shows the temperature profile for a short CPAIMD simulation of bulk silicon together with the kinetic energy profile from the fictitious orbital dynamics.The figure demonstrates that the orbital dynamics is essentially a "slave" to the nuclear dynamics, which shows that the electronic configuration closely follows that dynamics of the nuclei in the spirit of the Born-Oppenheimer approximation.

Plane wave basis sets and surface boundary conditions
In AIMD calculations, the most commonly employed boundary conditions are periodic boundary conditions, in which the system is replicated infinitely in all three spatial directions.This is clearly a natural choice for solids and is particularly convenient for liquids.In an infinite periodic system, the KS orbitals become Bloch functions of the form where k is a vector in the first Brioullin zone and u ik (r) is a periodic function.A natural basis set for expanding a periodic function is the Fourier or plane wave basis set, in which u ik (r) is expanded according to where V is the volume of the cell, g = 2πh −1 ĝ is a reciprocal lattice vector, h is the cell matrix, whose columns are the cell vectors (V = det(h)), ĝ is a vector of integers, and {c k i,g } are the expansion coefficients.An advantage of plane waves is that the sums needed to go back and forth between reciprocal space and real space can be performed efficiently using fast Fourier transforms (FFTs).In general, the properties of a periodic system are only correctly described if a sufficient number of k-vectors are sampled from the Brioullin zone.However, for the applications we will consider, we are able to choose sufficiently large system sizes that we can restrict our k-point sampling to the single point, k =( 0, 0, 0), known as the Γ-point.At the Γ-point, the plane wave expansion reduces to At the Γ-point, the orbitals can always be chosen to be real functions.Therefore, the plane-wave expansion coefficients satisfy the property that c * i,g = c i,−g ,w h i c hr e q u i r e s keeping only half of the full set of plane-wave expansion coefficients.In actual applications, plane waves up to a given cutoff |g| 2 /2 < E cut are retained.Similarly, the density n(r) given by Eq. ( 4) can also be expanded in a plane wave basis: However, since n(r) is obtained as a square of the KS orbitals, the cutoff needed for this expansion is 4E cut for consistency with the orbital expansion.
At first glance, it might seem that plane waves are ill-suited to treat surfaces because of their two-dimensional periodicity.However, in a series of papers (Minary et al., 2004;2002; 236 Silicon Carbide -Materials, Processing and Applications in Electronic Devices www.intechopen.com Tuckerman & Martyna, 1999), Martyna, Tuckerman, and coworkers showed that clusters (systems with no periodicity), wires (systems with one periodic dimension), and surfaces (systems with two periodic dimensions) could all be treated using a plane-wave basis within a single unified formalism.Let n(r) be a particle density with a Fourier expansion given by Eq. ( 19), and let φ(r − r ′ ) denote an interaction potential.In a fully periodic system, the energy of a system described by n(r) and φ(r − r ′ ) is given by where φg is the Fourier transform of the potential.For systems with fewer than three periodic dimensions, the idea is to replace Eq. ( 20) with its first-image approximation where φg denotes a Fourier expansion coefficient of the potential in the non-periodic dimensions and a Fourier transform along the periodic dimensions.For clusters, φg is given by φg = for wires, it becomes φg = and for surfaces, we obtain The error in the first-image approximation drops off as a function of the volume, area, or length in the non-periodic directions, as analyzed in Minary et al. (2004;2002); Tuckerman & Martyna (1999).
In order to have an expression that is easily computed within the plane wave description, consider two functions φ long (r) and φ short (r), which are assumed to be the long and short range contributions to the total potential, i.e.
We require that φ short (r) vanish exponentially quickly at large distances from the center of the parallelepiped and that φ long (r) contain the long range dependence of the full potential, φ(r).
with exponentially small error, ǫ(g), provided the range of φ short (r) is small compared size of the parallelepiped.In order to ensure that Eq. ( 26) is satisfied, a convergence parameter, α,is introduced which can be used to adjust the range of φ short (r) such that ǫ(g) ∼ 0andtheerror , ǫ(g), will be neglected in the following.
(28) Thus, Eq. ( 28) becomes leads to The new function appearing in the average potential energy, Eq. ( 29), is the difference between the Fourier series and Fourier transform form of the long range part of the potential energy and will be referred to as the screening function because it is constructed to "screen" the interaction of the system with an infinite array of periodic images.The specific case of the Coulomb potential, φ(r)=1/r, can be separated into short and long range components via where the first term is long range.The parameter α determines the specific ranges of these terms.The screening function for the cluster case is easily computed by introducing an FFT grid and performing the integration numerically (Tuckerman & Martyna, 1999).For the wire (Minary et al., 2002) and surface (Minary et al., 2004) When a plane wave basis set is employed, the external energy is made somewhat complicated by the fact that very large basis sets are needed to treat the rapid spatial fluctuations of core electrons.Therefore, core electrons are often replaced by atomic pseudopotentials or augmented plane wave techniques.Here, we shall discuss the former.In the atomic pseudopotential scheme, the nucleus plus the core electrons are treated in a frozen core type approximation as an "ion" carrying only the valence charge.In order to make this approximation, the valence orbitals, which, in principle must be orthogonal to the core orbitals, must see a different pseudopotential for each angular momentum component in the core, which means that the pseudopotential must generally be nonlocal.In order to see this, we consider a potential operator of the form where r is the distance from the ion, and |lm lm| is a projection operator onto each angular momentum component.In order to truncate the infinite sum over l in Eq. ( 32), we assume that for some l ≥ l, v l (r)=v l (r) and add and subtract the function v l (r) in Eq. ( 32): where the second line follows from the fact that the sum of the projection operators is unity, Δv l (r)=v l (r) − v l (r), and the sum in the third line is truncated before Δv l (r)=0.The complete pseudopotential operator is where v loc (r) ≡ v l (r) is known as the local part of the pseudopotential (having no projection operator attached to it).Now, the external energy, being derived from the ground-state expectation value of a one-body operator, is given by The first (local) term gives simply a local energy of the form which can be evaluated in reciprocal space as where Ṽloc (g) is the Fourier transform of the local potential.Note that at g =( 0, 0, 0),o n l y the nonsingular part of ṽloc (g) contributes.In the evaluation of the local term, it is often convenient to add and subtract a long-range term of the form Z I erf(α I r)/r,w h e r ee r f (x) is the error function, each ion in order to obtain the nonsingular part explicitly and a residual short-range function vloc (|r

Electron localization methods
An important feature of the KS energy functional is the fact that the total energy E[{ψ}, R] is invariant with respect to a unitary transformation within space of occupied orbitals.That is, if we introduce a new set of orbitals ψ ′ i (r) related to the ψ i (r) by where We say that the energy is invariant with respect to the group SU(N s ), i.e., the group of all N s × N s unitary matrices with unit determinant.This invariance is a type of gauge invariance, specifically that in the occupied orbital subspace.The fictitious orbital dynamics of the AIMD scheme as written in Eqs. ( 15) does not preserve any particular unitary representation or gauge of the orbitals but allows the orbitals to mix arbitrarily according to Eq. ( 38).This mixing happens intrinsically as part of the dynamics rather than by explicit application of the unitary transformation.
Although this arbitrariness has no effect on the nuclear dynamics, it is often desirable for the orbitals to be in a particular unitary representation.For example, we might wish to have the true Kohn-Sham orbitals at each step in an AIMD simulation in order to calculate the Kohn-Sham eigenvalues and generate the corresponding density of states from a histogram of these eigenvalues.This would require choosing U ij to be the unitary transformation that diagonalizes the matrix of Lagrange multipliers in Eq. ( 6).Another important representation is that in which the orbitals are maximally localized in real space.In this representation, the orbitals are closest to the classic "textbook" molecular orbital picture.
In order to obtain the unitary transformation U ij that generates maximally localized orbitals, we seek a functional that measures the total spatial spread of the orbitals.One possibility for this functional is simply to use the variance of the position operator r with respect to each orbital and sum these variances: The procedure for obtaining the maximally localized orbitals is to introduce the transformation in Eq. ( 38) into Eq.( 39) and then to minimize the spread functional with respect to The minimization must be carried out subject to the constraint the U ij be an element of SU(N s ).This constraint condition can be eliminated if we choose U to have the form U = exp(iA), where A is an N s × N s Hermitian matrix, and performing the minimization of Ω with respect to A.
A little reflection reveals that the spread functional in Eq. ( 39) is actually not suitable for periodic systems.The reason for this is that the position operator r lacks the translational invariance of the underlying periodic supercell.A generalization of the spread functional that does not suffer from this deficiency is (Berghold et al., 2000;Resta & Sorella, 1999) where σ and L denote the typical spatial extent of a localized orbital and box length, respectively, and Here G I = 2π(h −1 ) T ĝI ,w h e r e ĝI =( l I , m I , n I ) is the Ith Miller index and ω I is a weight having dimensions of (length) 2 .The function f (|z| 2 ) is often taken to be 1 −|z| 2 ,a l t h o u g h several choices are possible.The orbitals that result from minimizing Eq. ( 41) are known as Wannier orbitals |w i .I fz I,ii is evaluated with respect to these orbitals, then the orbital centers, known as Wannier centers, can be computed according to Wannier orbitals and their centers are useful in analyzing chemically reactive systems and will be employed in the present surface chemistry studies.Like the KS energy, the fictitious CP dynamics is invariant with respect to gauge transformations of the form given in Eq. (38).They are not, however, invariant under time-dependent unitary transformations of the form and consequently, the orbital gauge changes at each step of an AIMD simulation.If, however, we impose the requirement of invariance under Eq.( 44) on the CP dynamics, then not only would we obtain a gauge-invariant version of the AIMD algorithm, but we could also then fix a particular orbital gauge and have this gauge be preserved under the CP evolution.Using techniques for gauge field theory, it is possible to devise such a AIMD algorithm (Thomas et al., 2004).Introducing orbital momenta |π i conjugate to the orbital degrees of freedom, the gauge-invariant AIMD equations of motion have the basic structure where Here, the terms involving the matrix B ij (t) are gauge-fixing terms that preserve a desired orbital gauge.If we choose the unitary transformation U ij (t) to be the matrix that satisfies Eq. ( 40), then Eqs.(45) will propagate maximally localized orbitals (Iftimie et al., 2004).As was shown in Iftimie et al. (2004); Thomas et al. (2004), it is possible to evaluate the gauge-fixing terms in a way that does not require explicit minimization of the spread functional (Sharma et al., 2003).In this way, if the orbitals are initially localized, they remain localized throughout the trajectory.While the Wannier orbitals and Wannier centers are useful concepts, it is also useful to have a measure of electron localization that does not depend on a specific orbital representation, as the latter does have some arbitrariness associated with it.An alternative measure of electron localization that involves only the electron density n(r) and the so-called kinetic energy density was introduced by Becke and Edgecombe (1990).Defining the ratio χ(r)=D(r)/D 0 (r), where the function f (r)=1/(1 + χ 2 (r)) can be shown to lie in the interval f (r) ∈ [0, 1],where f (r)= 1 corresponds to perfect localization, and f (r)=1/2 corresponds to a gas-like localization.
The function f (r) is known as the electron localization function or ELF.In the studies to be presented below, we will make use both of the ELF and the Wannier orbitals and centers to quantify electron localization.

Reactions on the 3C-SiC(001)-3×2 surface
Silicon-carbide (SiC) and its associated reactions with a conjugated diene is an interesting surface to study and to compare to the pure silicon surface.In previous work (Hayes & Tuckerman, 2007;Minary & Tuckerman, 2004;2005), we have shown that when a conjugated diene reacts with the Si(100)-2×1 surface, a relatively broad distribution of products results, in agreement with experiment (Teague & Boland, 2003;2004), because the surface dimers are relatively closely spaced.Because of this, creating ordered organic layers on this surface using conjugated dienes seems unlikely unless some method can be found to enhance the population of one of the adducts, rendering the remaining adducts negligible.SiC exhibits a number of complicated surface reconstructions depending on the surface orientation and growth conditions.Some of these reconstructions offer the intriguing possibility of restricting the product distribution due to the fact that carbon-carbon or silicon-silicon dimer spacings are considerably larger.SiC is often the material of choice for electronic and sensor applications under extreme conditions (Capano & Trew, 1997;Mélinon et al., 2007;Starke, 2004) or subject to biocompatibility constraints (Stutzmann et al., 2006).Although most reconstructions are still being debated both experimentally and theoretically (Pollmann & Krüger, 2004;Soukiassian & Enriquez, 2004), there is widespread agreement on the structure of the 3C-SiC(001)-3×2 surface (D'angelo et al., 2003;Tejeda et al., 2004)(see Fig. 1), which will be studied in this section.SiC(001) shares the same zinc blend structure as pure Si( 001), but with alternating layers of Si and C. The top three layers are Si, the bottom in bulk-like positions and the top decomposed into an open 2/3 + 1/3 adlayer structure.Si atoms in the bottom two-thirds layers are 4-fold coordinated dimers while those Si atoms in the top one-third are asymmetric tilted dimers with dangling bonds.Given the Si-rich surface environment and presence of asymmetric surface dimers, one might expect much of the same Si-based chemistry to occur with two significant differences: (1) altered reactivity due to the surface strain (the SiC lattice constant is ∼ 20% smaller than Si) and (2) suppression of interdimer adducts due to the larger dimer spacing compared to Si (∼60% along a dimer row, ∼20% across dimer rows).Previous theoretical studies used either static (0 K) DFT calculations of hydrogen (Chang et al., 2005;Di Felice et al., 2005;Peng et al., 2007a;2005;2007b), a carbon nanotube (de Brito Mota & de Castilho, 2006), or ethylene/acetylene (Wieferink et al., 2006;2007) adsorbed on SiC(001)-3×2 or employed molecular dynamics of water (Cicero et al., 2004) or small molecules of the CH 3 -X family (Cicero & Catellani, 2005) on the less thermodynamically stable SiC(001)-2×1 surface.Here, we consider cycloaddition reactions on the SiC-3×2 surface that include dynamic and thermal effects.A primary goal for considering this surface is to determine whether 3C-SiC(001)-3×2 is a promising candidate for creating ordered semiconducting-organic interfaces via cycloaddition reactions.In the study Hayes and Tuckerman (2008), the KS orbitals were expanded in a plane-wave basis set up to a kinetic energy cut-off of 40 Ry.As in the 1,3-CHD studies described above, exchange and correlation are treated with the spin restricted form of the PBE functional (Perdew et al., 1996), and core electrons were replaced by Troullier-Martins pseudopotentials (Troullier & Martins, 1991) with S, P, and D treated as local for H, C, and Si, respectively.The resulting SiC theoretical lattice constant, 4.39 Å, agrees well with the experimental value of 4.36 Å (Tejeda et al., 2004).The full system is shown in Fig. 1.The 3×2 unit cell is doubled in both directions to include four surface dimers to allow the possibility of all interdimer adducts.Again, the resulting large surface area, (18.6 Å x 12.4 Å), allows the Γ-point approximation to be used in lieu of explicit k-point sampling.Two bulk layers of Si and C, terminated by H on the bottom surface, provide a reconstructed (1/3 + 2/3) Si surface in reasonable agreement with experiment (see below).The final system has 182 atoms [24 atoms/layer * (1 Si adlayer + 4 atomic layers) + 2*24 terminating H].The simulation cell employed lengths of 18.6 Å and 12.4 Å along the periodic directions and 31.2Å along the nonperiodic z direction.Both the CHD and SiC(001) surface were equilibrated separately under NVT conditions using Nosé-Hoover chain thermostats (Martyna et al., 1992) at 300 K with a timestep of 0.1 fs for 1 ps and 3 ps, respectively.When the equilibrated CHD was allowed to react with the equilibrated surface, the time step was reduced to 0.05 fs in order to ensure adiabaticity.The CHD was placed 3 Å above the surface, as defined by the lowest point on the CHD and the highest point on the surface.Each of twelve trajectories was initiated from the same CHD and SiC structures but with the CHD placed at a different orientations and/or locations over the surface.The subsequent initialization procedure was identical to the CHD-Si(100) system: First the system was annealed from 0 K to 300 K in the NVE ensemble.Following this, it was equilibrated with Nosé-Hoover chain thermostats for 1 ps at 300K under NVT conditions, keeping the center of mass of the CHD fixed.Finally, the CHD center of mass constraint was removed and the system was allowed to evolve under the NVE ensemble until an adduct formed or 20 ps elapsed.The reactions that occur on this surface all take place on or in the vicinity of a single surface Si-Si dimer.However, as Fig. 2 shows, there is not one but rather four adducts that are observed to form.Adduct labels from the Si + CHD study are used for consistency.As postulated, the widely spaced dimers successfully suppressed the interdimer adducts that formed on the Si(100)-2×1 surface (Hayes & Tuckerman, 2007).From the twelve trajectories, three formed the [4+2] Diels-Alder type intradimer adduct (A), one produced the [2+2] intradimer adduct (D), five exhibited hydrogen abstraction (H), and one resulted in a novel [4+2] subdimer adduct between Si in d 1 and d 2 (G) (see Fig. 1).The remaining trajectories only formed 1 C-Si bond within 20 ps.Although the statistics are limited, these results suggest that H abstraction is favorable, consistent with the high reactivity of atomic H observed in experimental studies on this system (Amy & Chabal, 2003;Derycke et al., 2003).What is somewhat more troublesome, from the point of view of creating well-ordered  C, and H are represented by yellow, blue, and silver, respectively.The remaining C=C bond(s) is highlighted in green.The larger spacing between dimers suppresses interdimer adducts.However, adduct (G) destroys the surface, rendering this system inappropriate for applications requiring well-defined organic-semiconducting interfaces.
organic-semiconducting interfaces is the presence of the subdimer adduct G.All the surface bonds directly connected to the adduct slightly expand to 2.42-2.47Å, with the exception of one bond to a Si in the third layer, (highlighted in red in Fig. 2G) which disappears entirely.The energetic gain of the additional strong C-Si bond outweighs the loss of a strained Si-Si bond.The end effect is the destruction of the perfect surface and the creation of an unsaturated Si in the bulk.One adduct is noticeably missing: the [2+2] subdimer adduct.At several points during the simulation this adduct was poised to form but quickly left the vicinity.Most likely, the strain caused by the four-member ring combined with the two energetically less stable unsaturated Si prevented this adduct from forming, even though the [2+2] intradimer and [4+2] subdimer adducts are stable.In Fig. 3, we show the carbon-carbon and CHD-Si distances as functions of time for the different adducts observed.This figure reveals that the mechanism of the reactions proceeds in a manner very similar to that of CHD and 1,3-butadiene on the Si(100)-2×1surface:Itisan asymmetric, nonconcerted mechanism that involves a carbocation intermediate.What differs from Si(100) is the time elapsed before the first bond forms and the intermediate lifetime.On the Si(100)-2×1 surface the CHD always found an available "down" Si to form the first bond within less than 10 ps or 40 Å of wandering over the surface.On the SiC(001)-3×2surfacethe exploration process sometimes required up to 20 ps and over 100 Å.While the exact numbers are only qualitative, the trend is significant.The Si(100) dimers are more tilted on average, and hence expected to be slightly more reactive.However, the dominant contribution is likely the density of tilted dimers: Si has 0.033 dimers/Å 2 , but SiC only has 0.017 dimers/Å 2 .Regardless of whether dimer flipping occurs, it is simply more difficult to find a dimer on the SiC surface.An important consideration in cycloaddition reactions such as those studied here is the possibility of their occurring through a radical mechanism.Multi-reference self consistent field cluster calculations of the SiC(001)-2x1 surface suggest that the topmost dimer exhibits significant diradical character (Tamura & Gordon, 2003), and since DFT is a single-reference method, multi-reference contributions are generally not included.
However, cluster methods may bias the results by unphysically truncating the system instead of treating the full periodicity.For instance, cluster methods predict that Si(100)-2×1d i m e r sa r e symmetric (Olson & Gordon, 2006), contrary to experimental evidence (Mizuno et al., 2004;Over et al., 1997), while periodic DFT correctly captures the dimer tilt (Hayes & Tuckerman, 2007).In order to estimate the importance of diradical mechanisms and surface crossing, a series of single point energy calculations at regular intervals during four representative trajectories are plotted in Fig. 4. Three electronic configurations are considered: singlet spin restricted (SR) where the up and down spin are identical (black down triangles), singlet spin unrestricted (SU) where the up and down spin can vary spatially (red up triangles), and triplet SU (green squares).In all cases, the triplet configuration is unfavorable.However, at two places in the transition state (Adduct A in Fig. 4a at 3000 fs and Adduct G in Fig. 4d at 8500 fs) the single SU is slightly favored.Thus, multi-reference methods, which can account for surface crossing, may yield alternative reaction mechanisms.subdimer trajectory.Energies are relative to the value at t=0 in (a).The insets show the difference between the spin restricted and unrestricted singlet energies (blue line) and spin restricted and triplet unrestricted energies (purple).The triplet configuration is always highest in energy.In (a) at 3000 fs and (d) at 8500 fs the singlet transition state configuration is favorable by 2.3 and 0.5 kcal/mol, respectively.Thus, a radical mechanism may also occur in this system.which we take to be the Diels-Alder-type [4+2] intradimer product.In the present case, the expect some of the barriers to product formation to be sufficiently high that specialized free energy sampling techniques are needed.Here, we employ the so-called blue moon ensemble approach (Carter et al., 1989;Sprik & Ciccotti, 1998) combined with thermodynamic integration.In order to define such a free energy profile, we first need to specify a reaction coordinate capable of following the progress of the reaction.For this purpose, we choose a coordinate ξ of the form where C m 1 and C m 4 are the carbons in the 1 and 4 positions in the organic molecule, and C s a and C s b are the two surface carbon atoms with which covalent bonds will form with C m 1 and C m 4 .Over the course of the reactions considered, ξ decreases from approximately 4 Å to a value less than 1.5 Å.In the aforementioned blue moon ensemble approach (Carter et al., 1989;Sprik & Ciccotti, 1998), the coordinate ξ is constrained at a set of equally spaced points between the two endpoints.At each constrained value, an AIMD simulation is performed over which we compute a conditional average ∂H/∂ξ cond ,w h e r eH is the nuclear Hamiltonian.Finally, the full free energy profile is reconstructed via thermodynamic integration: An example of such a free energy profile for the [4+2] cycloaddition reaction of 1,3-butadiene with a single silicon surface dimer on the Si(100)-2×1s u r f a c ei ss h o w ni n Fig. 6 (Minary & Tuckerman, 2004).We show this profile as an example in order that direct comparison can be made between reactions on this surface and those on the SiC(100)2×2 surface.The profile in Fig. 6 shows an initial barrier at ξ =3.In the present study of 1,3-cyclohexadiene with the SiC(100)-2×2 surface, the KS orbitals were expanded in a plane-wave basis set up to a kinetic energy cut-off of 60 Ry.As in the 1,3-CHD studies described above, exchange and correlation are treated with the spin restricted form of the PBE functional (Perdew et al., 1996), and core electrons were replaced by Troullier-Martins pseudopotentials (Troullier & Martins, 1991) with S, P, and D treated as local for H, C, and Si, respectively.The periodic slab contains 128 atoms arranged in 6 layers (including a bottom passivating hydrogen layer).Proper treatment of surface boundary conditions allowed for a simulation cell with dimensions L x = 17.56Å,L y =8.78 Å, and L z =31 Å along the nonperiodic dimension.The surface contains 8 C≡C dimers.This setup is capable of reproducing the experimentally observed dimer buckling (Derycke et al., 2000) that static ab initio calculations using cluster models are unable to describe (Bermudez, 2003).In Fig. 7, we show the free energy profile for the [4+2] cycloaddition reaction of 1,3-cyclohexadiene with one of the C≡C surface dimers.The free energy profile is calculated by dividing the ξ interval ξ ∈ [1.59, 3.69] into 15 equally spaced intervals, and each constrained simulation was equilibrated for 1.0 ps followed by 3.0 ps of averaging using a time step of 0.025 fs.All calculations are carried out in the NVT ensemble at 300 K using Nosé-Hoover chain thermostats (Martyna et al., 1992) In contrast to the free energy profile of Fig. 6, the profile in Fig. 7 shows no evidence of a stable intermediate.Rather, apart from an initial barrier of approximately 8 kcal/mol, the free energy is strictly downhill.The reaction is thermodynamically favored by approximately 48 kcal/mol.The suggestion from Fig. 7 is that the reaction is symmetric and concerted in contrast to the reactions on the other surfaces we have considered thus far.Fig. 7 shows snapshots of the molecule and the surface atoms Further evidence for the concerted asymmetric nature of the reaction is provided in Fig. 8, which shows the average carbon-carbon lengths computed over the constrained trajectories at each point of the free energy profile.It can be seen by the fact that one CC bond forms before the other that there is a slight tendency for an asymmetric reaction, despite its being concerted.
In order to demonstrate that the [4+2] Diels-Alder type cycloaddition product is highly favored over other reaction products on this surface, we show one additional example of a free energy profile, specifically, that for the formation of a [2+2] cycloaddition reaction with a single surface C≡C dimer.This profile is shown in Fig. 9.In contrast to the [4+2] Diels-Alder type adduct, the barrier to formation of this adduct is roughly 27 kcal/mol (compared to 8 kcal/mol for the Diels-Alder product).Thus, although the [2+2] reaction  is thermodynamically favorable, this barrier is sufficiently high that we would expect this particular reaction channel to be highly suppressed compared to one with a 19 kcal/mol lower barrier.The free energy profile, together with the snapshots taken along the reaction path, also suggests that this reaction occurs via an asymmetric, concerted mechanism, as was found for the [4+2] Diels-Alder type product.Although we have not shown them here, we have computed free energy profiles for a variety of other adducts, including interdimer and sublayer adducts, and in all cases, free energy barriers exceeding 20 kcal/mol (or 40 kcal/mol in the case of the sublayer adduct) were obtained.These results strongly suggest that the product distribution this surface would, for all intents and purposes, restricted to the single [4+2] Diels-Alder type product, implying that this surface might be a candidate for creating an ordered organic/semiconductor interface.

237
Creation of Ordered Layers on Semiconductor Surfaces: An ab Initio Molecular Dynamics Study of the SiC(001)-3x2 and SiC(100)-c(2x2) Surfaces www.intechopen.comWith these two requirements, it is possible to write φshort

Fig. 3 .
Fig. 3. Relevant bond lengths () vs time (fs) during product formation for four representative adducts.The top row displays the C-C bonds lengths (moving average over 25 fs) while the bottom row plots the first and second CHD -SiC surface bond.The color-coded inset identifies the bond being plotted for each adduct type.Change in the C-C bond length closely correlates with surface-adduct bond formation.Intermediate lifetimes over all trajectories range from 0.05 -18 + ps.

Fig. 4 .
Fig. 4. Spin restricted (black down triangles), singlet spin-unrestricted (red up triangles), and triplet unrestricted (green squares) energies at regular intervals during a representative (a) [4+2] intradimer adduct [A], (b) [2+2] intradimer adduct [D], (c) H-abstraction, and (d) [4+2]subdimer trajectory.Energies are relative to the value at t=0 in (a).The insets show the difference between the spin restricted and unrestricted singlet energies (blue line) and spin restricted and triplet unrestricted energies (purple).The triplet configuration is always highest in energy.In (a) at 3000 fs and (d) at 8500 fs the singlet transition state configuration is favorable by 2.3 and 0.5 kcal/mol, respectively.Thus, a radical mechanism may also occur in this system.

Fig. 6 .
Fig. 6.Free energy along the reaction pathway leading to a Diels-Alder [4+2] adduct.Blue and red triangles indicate the product (EQ) and intermediate states (IS), respectively.Inset shows the buckling angle (α) distribution of the Si dimer for both the IS (red) and the EQ configurations (blue).The snapshots include configurations representing the IS and EQ geometries.Blue, green, and white spheres denote Si, C, and H atoms, respectively, and gray spheres indicate the location of Wannier centers.Red spheres locate positively charged atoms.The purple surface is a 0.95 electron localization function (ELF) isosurface.

Fig. 7 .
Fig. 7. Free energy profile for the formation of the [4+2] Diels-Alder-like adduct between 1,3-cyclohexadiene a C≡C dimer on the SiC-2×2 surface.Blue, white and yellow spheres represent C, H, and Si, respectively.Red spheres are the centers of maximally localized Wannier functions.withwhich it interacts at various points along the free energy profile.In these snapshots, red spheres represent the centers of maximally localized Wannier functions.These provide a visual picture of where new covalent bonds are forming as the reaction coordinate ξ is decreased.By following these, we clearly see that one CC bond forms before the other, demonstrating the asymmetry of the reaction, which is a result of the buckling of the surface dimers.The buckling gives rise to a charge asymmetry in the C≡Cs u r f a c ed i m e r ,a n da s a result, the first step in the reaction is a nucleophilic attack of one of the C=Cb o n d si n the cyclohexadiene on the positively charged carbon in the surface dimer, this carbon being the lower of the two.Once this first CC bond forms, the second CC bond follows after a change of approximately 0.3 Å in the reaction coordinate with no stable intermediate along the way toward the final [4+2] cycloaddition product.In addition, the Wannier centers show the conversion of the triple bond on the surface to a double bond in the final product state.Further evidence for the concerted asymmetric nature of the reaction is provided in Fig.8, which shows the average carbon-carbon lengths computed over the constrained trajectories at each point of the free energy profile.It can be seen by the fact that one CC bond forms before the other that there is a slight tendency for an asymmetric reaction, despite its being concerted.In order to demonstrate that the [4+2] Diels-Alder type cycloaddition product is highly favored over other reaction products on this surface, we show one additional example of a free energy profile, specifically, that for the formation of a [2+2] cycloaddition reaction with a single surface C≡C dimer.This profile is shown in Fig.9.In contrast to the [4+2] Diels-Alder type adduct, the barrier to formation of this adduct is roughly 27 kcal/mol (compared to 8 kcal/mol for the Diels-Alder product).Thus, although the [2+2] reaction

Fig. 9 .
Fig. 9. Free energy profile for the formation of the [2+2] adduct between 1,3-cyclohexadiene a C≡C dimer on the SiC-2×2 surface.Blue, white and yellow spheres represent C, H, and Si, respectively.Red spheres are the centers of maximally localized Wannier functions.

2
(Hayes & Tuckerman, 2007;/mol.As ξ decreases, a shallow minimum/plateau is seen at ξ =2.75 Å, and such a minimum indicates a stable intermediate.This intermediate was identified as a carbocation in which one of the C-Si bonds had formed prior to the second bond formation(Hayes & Tuckerman, 2007;