Contents In this project, DFT calculations have been

Contents Project 1. DFT – Structural stability and heat of formation 1. Introduction 32. Convergence Criteria 43. Results and Discussion 53.1 Structural stability 53.2 Density of states of Si (fcc) 83.3 Heat of formation 93.3.1 Formation energy of AlSi 93.3.2 Formation energy of TiAl 94. Multiscale approach 105. Summary 10  Project 3. KMC – simple model of a catalytic reaction 1. Introduction 112. Results and Discussion 122.1 Diffusion of a single atom 122.2 Reaction rate 142.2.1 Concentration dependence 142.2.2 Temperature dependence 152.2.3 Formulation of reaction rate using KMC model 193. Future Outlook 204. Summary 20 

Project 1. DFT – Structural stability and heat of formation1. IntroductionHohenberg and Kohn laid the theoretical foundations of density functional theory (DFT) in 1964, proving that the electron density uniquely determines the total energy of the electronic ground state in a many-body system. Assuming the Born-Oppenheimer approximation, this allows a set of equations to describe the behaviour of non-interacting electrons in an effective local potential, known as the Kohn-Sham potential:The Kohn-Sham potential ?s(r) is given by the following: where r is the electron’s position in space, vext is external potential, n(r) is the electron density, vXC is the exchange-correlation term, and EXC is the exchange-correlation function. The exact form of the exchange-correlation function is unknown. Therefore, the reliability of DFT calculations depends on the approximations for the exchange-correlation energy functional. The difficulty in plan-wave expansion of the Kohn-Sham equations arises from the rapid oscillations in wave function inside the core region. By using pseudopotentials, the eigenvalues of the valence states remain unchanged, whereas the wave function has been smoothed and can be efficiently solved through self-consistent calculations. specified convergence criteria are required for total energy calculation over points in k-space. Ground-state properties such as charge density, total energy, elastic constants etc. obtained are generally reliable estimates. DFT provides us with a relatively efficient and unbiased tool to compute ground-state properties and electronic structure of a system by directly solving approximate versions of the Schrödinger equation. In this project, DFT calculations have been performed on different structures of Silicon, Aluminium, Titanium and their alloys. All calculations were carried out by using the Cambridge Sequential Total Energy Package (CASTEP) with ultrasoft pseudo-potentials and a plane-wave basis-set. The exchange and correlation energies were calculated within the local-density approximation (LDA). The energies obtained are expressed in the same arbitrary scale, which allow the structural stability and heat of formation to be determined and compared. 2. Convergence CriteriaThe convergence tests performed for total energy of silicon with respect to (w.r.t.) k-point density and w.r.t. cutoff energy, illustrated in Fig. 1, show that the results were converged to within 1 mÅ for k value greater than 4, and to within 4 mÅ for Ecutoff larger than 200 eV. Therefore k-point of 10 and Ecutoff of 300eV were used for Si calculation to produce reliable results while saving the computational time.  Figure 1. Convergence of lattice constant of Si (diamond) versus Ecutoff and k-point density Metals are much more sensitive to k-point sampling, because for some bands, a part of the Brillouin zone is filled and the remainder is empty, separated by a Fermi surface. Titanium, a transition metal, has partially filled d band and thus needs larger number of k points for computation. The number of k-points directly determines how accurately the Fermi surface(s) can be resolved. Figure 2. Convergence of lattice constant of Ti (bcc) versus Ecutoff and k-point density Convergence tests for Ti in body-centred cubic structure were conducted and the results are plotted in Fig. 2. The lattice constants were converged to within 11 mÅ for k > 6 and to within 2 mÅ for Ecutoff > 350 eV. Thus, k of 10 and Ecutoff of 450 eV were chosen for Ti calculation.    3. Results and Discussion3.1 Structural stabilityTotal energies of four types of cubic structures, namely diamond, body-centered cubic (bcc), face-centered cubic (fcc), simple cubic (sc), were calculated by means of CASTEP in an arbitrary scale and compared to find the most thermodynamically stable phase for Si, Al and Ti. For compounds AlSi and TiAl, the cubic B1 (NaCl), B2 (CsCl), B3 (ZnS), and tetragonal L10 (CuAu) structure types were used to determine the relative stability. The results are portrayed in the graphs below (Fig. 3-7) and summarised in the following table.    Figure 3. The calculated total energies per atom of silicon as a function of lattice constant    Figure 4. The calculated total energies per atom of aluminium as a function of lattice parameter     Figure 5. The calculated total energies per atom of titanium as a function of lattice parameter  Figure 6. The calculated total energies per atom of AlSi as a function of lattice parameter     Figure 7. The calculated total energies per atom of TiAl as a function of lattice parameter   Element/CompoundStable phaseLowest Etot/eV/atomCorresponding Lattice constant/ÅSidiamond-165.45.40Alfcc-76.43.98Tifcc-1600.84.00AlSiL10-118.93.88TiAlL10-839.13.94 Table 1. Calculated lowest total energies, lattice parameters and the most stable structures of listed elements and compounds 3.2 Density of states of Si (fcc)As demonstrated in Fig. 8, Si in diamond structure exhibits insulating behavior with a band gap, where density of states (DOS) reaches zero for a range of energies around Fermi energy. It should be noted that DFT within the LDA approximation is known to underestimate band gaps and as a consequence, the actual band gap would be larger than what is shown here.  Figure 8. DOS versus energy of states of Si in diamond structure However, silicon in fcc structure shows metallic behaviour with no band gap present. As having similar electronic configurations with valence electrons in 3s and 3p orbitals, the pseudopotentials used for Si (fcc) and Al, which accounts for the repulsion from inner ion core would be close to each other. Thus, DOS plots of Si (fcc) and Al calculated under this approximation are in a similar pattern.  Figure 9. DOS versus energy of states of (a) Si in fcc structure and (b) Al The sudden drop in DOS is due to the pseudopotential gap between bands and the upper limit in the energy scale is determined by the number of bands generated from pseudowavefunctions. In comparison to Fig. 9 (b), the Fermi level shifts to a higher value in Fig. 9 (a) since Si has one more valence electron to fill in the states than Al.  3.3 Heat of formation3.3.1 Formation energy of AlSiThe heat of formation per atom of AlSi can be determined from where s stands for the energetically most stable phase. Therefore,The formation of a molecule of AlSi would require two atoms, one Al atom and one Si atom, which means that the relative formation energy of AlSi is +4.0 eV in total. It suggests that AlSi alloy would tend to decompose and rarely be observed in reality. 3.3.2 Formation energy of TiAlThe heat of formation per atom of TiAl can be determined from  Therefore,The formation energy of TiAl alloy is negative, which implies that Al can be dissolved into the lattice of Ti to form the compounds of TiAl. Thus, the alloy phase is expected to exist as it is more thermodynamically stable than individual Ti or Al atoms.  4. Multiscale approachDFT is effectively utilised as a proficient quantum mechanics method. Its applications range from studying the electronic structure, dynamics, and excited states to bulk modulus, dielectric constants and magnetic properties. However, picking the ‘right’ approximate exchange-correlation functional can be problematic, which can strongly affect the calculated results and might lead to errors.  For large systems, low scaling of DFT becomes a problem. There are two viable solutions to this: DFT calculation can be embedded into a framework considering both the local character of bonds and extended nature of electronic states in the bulk, or one can utilize a multiscale approach that is capable of combining speed and accuracy with the aim of capturing the majority of long-range effects. Nevertheless, DFT is a general tool for studying the many-body system at present, thanks to the high accuracy that can be achieved at relatively low computational cost.  5. SummaryIn this project, structural stability, density of states and heat of formation of given elements and compounds were investigated by using the density functional theory within LDA and ultrasoft pseudopotential approximations. It is found that diamond structure is most stable for silicon, while the fcc structure is most stable for Al and Ti. As for compounds, AlSi and TiAl would both adopt the L10 structure. Calculations of relative formation energies were carried out using the lowest total energies of the stable phases. A positive heat of formation suggests the structural instability of AlSi, while TiAl, with negative formation energy, is likely be found as a stable compound.        Project 3. KMC – simple model of a catalytic reaction1. Introduction1.1 The kinetic Monte Carlo methodThe kinetic Monte Carlo (KMC) method is a stochastic approach for modelling a sequence of thermally activated events such that the likelihood of each event happens is consistent with the rate of event. A typical KMC algorithm executes elementary events one after the other and generates a sequence of configurations and times to determine when the transitions between these configurations take place. This solves the master equations numerically to obtain accurate estimates: where ci(t) is the concentration of adatoms on site i, ?ji denotes the rate constant for adatom hops from i to j and the sum over j runs over all neighbouring configurations of lattice site i.  1.2 Reaction modelA simple model which describes a catalytic reaction between two different species A and B is used here. The catalyst surface is given by a two-dimensional integer lattice. The species A and B may diffuse from one adsorption site to another by overcoming energy barriers, given by ?EA and ?EB respectively. The atoms are assumed to react with each other instantly when they occupy the same adsorption site. Once the product A-B desorbs from the catalyst surface, the site will be re-occupied by incoming atoms.  Figure 10. Snapshot of diffusing adatoms of A and B (pink and orange spheres), and reaction product (blue sphere) on the surface of catalyst at kT = 1.0 This model assumes that the partial pressure can be modified to keep the total numbers of A and B constant. In practice, however, surface coverages are variables of a differential equation system to be solved for. Furthermore, all catalytic sites are considered to be equivalent, which is mostly untrue: catalytic surface exposes sites with various coordination numbers, onto which adsorbates bind with different strengths. And there is no information on how fast the desorption process of the product A-B is and the possibility of them kept fixed to the surface, which might cause reduction in number of available adsorption sites as reaction proceeds. These assumptions pose limitations on the validity of this model in mimicking realistic catalytic systems.  2. Results and Discussion2.1 Diffusion of a single atomThe diffusion of a single species was analysed by putting only one atom of A into a surface system (100×100 sites). Total number of KMC attempt steps conducted was set to be 107. The diffusion coefficient is defined as /?t, where is the mean square displacement from the initial position and ?t is the time between two measurements. Diffusion coefficients of species A were collected and the averaged results were compared. The relationship between diffusion coefficient of A and ?t is studied by changing the number of attempt steps per absorbed atom after which the current configuration is analysed (nmeasure). As demonstrated in Fig. 11, the diffusion coefficient is approximately constant over small ?t, but the deviation diverges as ?t increases. Figure 11. Diffusion coefficient versus log10(?t). Different nseeds were used to generate sequences of random numbers.The A atom would undergo a ‘random walk’ and move by hopping to neighbouring sites along trajectories generated by KMC algorithm. Therefore, the diffusion coefficients obtained from the mean square displacements should be roughly the same regardless of any changes in ?t, due to the equal probability of movement in all directions. When ?t is increased, the measurement is taken over a larger fraction of the total attempts and fewer measurements are made. The diffusion coefficient collected is only averaged over a small number of measurements, leading to large data scattering. The Arrhenius equation can be used to describe the diffusion process here:where v is the attempt frequency (pre-exponential factor) and Ea is the activation energy of diffusion. To get v and Ea, the diffusion coefficients for a range of temperatures were calculated and linear fit was done within the assumption that Ea and v are independent of temperature:The y-intercept would give the attempt frequency and the gradient of the linear fit equals to the activation energy. Simulation results and their linear regression are reported in Fig. 13.  Figure 13. ln(DA) versus 1/kT (in arbitrary unit)    In this case,  and . The value of activation energy, which is very close to the diffusion barrier of species A (?EA = 1.0), implies that surface diffusion of A atoms is the only process happened during the simulation.Therefore, diffusion coefficient of species A . 2.2 Reaction rate2.2.1 Concentration dependenceTwo series of simulations were carried out to investigate the relationship between reaction rate and surface concentration of reactants at a given temperature. In the first, concentration of B atoms was held constant at 25% to ensure that all A atoms react. In the second, the concentration of species B was set equal to that of A and the reaction rate was thus affected by two variables.   Figure 14. Concentration dependence of reaction rate on surface concentration (?EA = ?EB) Figure 15. Concentration dependence of reaction rate on surface concentration (?EA = 3?EB) As can be seen from Fig. 14 and Fig. 15, the reaction rate shows linear dependence on concentration when the other variables were fixed (despite some scattering at low concentration) and quadratic relationship when cA and cB were both varied. This can be explained by the rate law. In this model, A and B both have first order of reaction and hence the reaction rate can be expressed as:where k is the rate constant. For cB = 25%, For cA = cB, By comparing the equations of best fit, it is evident that the calculated results agree well with the rate law. 2.2.2 Temperature dependenceTo investigate the role of surface diffusion in the reaction, the diffusion barrier of species A was set to be 1.0 and 3.0 while keeping that of species B constant at 1.0 (in arbitrary unit). Simulations were carried out at surface coverages of cA = cB =5% and two Arrhenius plots were obtained. Figure 16. The logarithmic reaction rate vs the reciprocal of kT for ?EA = ?EB   Figure 17. ln(Reaction rate) vs 1/kT for ?EA = 3?EB  When the diffusion barriers of both species are equal, the gradient of fitted line is independent of temperature. By relating it to the Arrhenius equation, the activation energy can be determined:where A is the pre-exponential factor and Ea is the activation energy of reaction. In this case, the activation energy equals to 1.0013. It is apparent that the activation energy of reaction itself (0.0013) is negligible compared to the diffusion barrier (?EA = ?EB = 1.0). Hence, this type of reaction can be described as ‘diffusion-controlled’.  Figure 18. Illustrative diagram of the reaction event. The transitions of A and B atoms from one potential minimum to the next are equally likely as ?EA = ?EB = 1.0. When the diffusion of species A needs two times more energy than that of species B, the activation energy becomes dependent on temperature. The reaction rate of the following form is considered:  where aL is the pre-exponetial factor at low temperature, ?EL is the activation energy at low temperature, and aH is the pre-exponetial factor at high temperature ?EH is the activation energy at high temperature. At low temperature, kT ? 0,      The Taylor expansion of  =  (to the first order).So   At high temperature, ?0 The Taylor expansion of  =  (to the first order).Hence,  By solving the equations above, the pre-exponential factors and activation energies are evaluated and shown in Table 2.  ?EaLow T1.070.0149High T2.720.0076 Table 2. The activation energies and pre-exponential factors at temperatures below (Low T) and above (High T) kT = 0.74 Figure 19. Illustrative diagram of the reaction event. The probability of (a) B hopping to A is greater than that of (b) A hopping to B because of the difference in diffusion barriers. At low temperatures, the dominant process is B atoms jumping to neighbouring sites, which can be proven by comparing ?EL (1.07) and ?EB (1.0). At higher temperature, the diffusion of species A becomes significant due to higher kinetic energy the atoms possess. Thus, the activation energy at high temperatures is a weighted average of both diffusion barriers.   2.2.3 Coarse-grained KMC A formula can be extracted from the above analysis of effects of concentrations and temperatures on the catalytic reaction rate, which is a combination of Arrhenius and rate equations in chemistry.  where k is the rate constant, cA and cB are concentrations of species A and B, kL and kH are pre-exponential factors at low and high temperatures, ?EL and ?EH are the activation energies of reaction at low and high temperatures, which are related to diffusion barriers.  3. Future Outlook The KMC simulations that have been discussed focus on the interaction between reactants and a catalytic surface. Combination with other models would extend the KMC simulations to include  more processes dealing with transfer of matter and heat in real catalytic systems. In addition, DFT and molecular dynamics can be used to calculate exact values of attempt frequencies and energy barriers. This project can also be embedded into a continuum model, using finite element analysis to optimise the use of catalyst by controlling its operating conditions, such as working temperature, distribution and surface contact between catalyst and reactants.   4. SummaryThe KMC method has been a vital tool in advancing our understanding of catalytic systems. In this project, a surface catalysed reaction model of two species was developed and applied in KMC simulation. The results are in agreement with the Arrhenius equation and rate law in chemistry for both single-species diffusion and catalytic reaction. It is also shown that the diffusion process with lower energy barrier dominates in low-temperature reaction, whereas the one requires higher activation energy only becomes important upon heating. Moreover, a formula describing the combined effect of concentration and temperature on reaction rate has been suggested using  simulated results.