SCC-DFTB for Transition Metals
Tight Binding Parametrisation
By: Simon Colbin
Supervisor: Peter Broqvist
Abstract
In this report, Slater-Koster (SK) tables have been generated for the self-consistent charge density func- tional based tight binding method (SCC-DFTB) to accurately describe the electronic structure of tran- sition metals in group 10 and 11. The SK-tables were tuned by changing the compression radius, both for the potential and the wave-function, so that the resulting band structure calculated from bulk models using SCC-DFTB matched the band structure generated from the same model using density functional theory (DFT). The thesis work has resulted in 6 unique sets of parameters, one set for each of the con- sidered elements Ni, Cu, Pd, Pt, Ag, and Au. Successful tests of the transferability of the generated SK-tables was also performed. In particular, the d-band center of a 79 atom Ni nanoparticle was cal- culated using both the SCC-DFTB method and the corresponding DFT reference method showing an error less than 0.2 eV. In final, a demonstration of the potential use of the developed methodology is also presented, where the electronic structure of a Ni surface consisting of over 1000 atoms with several defects, such as steps and kinks, and vacancies of different sizes, were calculated using SCC-DFTB.
Contents
1 Introduction 1
2 Aim 3
3 Theory 4
3.1 Molecular orbitals . . . . 4
3.2 Band structure. . . . 5
3.3 Density of states. . . . 6
3.4 The d-band model . . . . 6
4 Method 8 4.1 DFT . . . . 8
4.2 SCC-DFTB . . . . 8
4.2.1 Optimisation . . . . 9
4.3 Calculating the d-band center . . . . 9
5 Results and Discussion 11 5.1 DFT . . . 11
5.2 The influence of the compression radius and electronic configuration on the band . . . . 11
5.3 Finding parameters using PSO . . . 16
5.4 Extending the parametrisation on other elements. . . 16
5.5 Testing transferability: d-band center for different Ni structures . . . 18
6 Conclusion 21
Abbreviations
ASE Atomic Simulation Environment DFT Density functional theory DFTB Density functional tight binding DOS Density of states
GGA-PBE Generalised gradient approximation Perdew-Burke-Ernzerhof
MO Molecular orbital
pDOS Projected density of states PSO Particle swarm optimization
SCC-DFTB Self-consistent charge density functional based tight binding method
SK Slater-Koster
VASP Vienna Ab initio Simulation Package
1 Introduction
Heterogeneous catalysis is extremely important in an industrial context, and has been a research topic for over a hundred years. An example that illustrates its importance is the Haber Bosch process where ammonia is created by catalysing a reaction between nitrogen and hydrogen on an iron surface. [1] In 1912, Paul Sabatier was awarded the Nobel prize for the use of catalysis in organic synthesis. He dis- covered that nickel could be used to improve the hydrogenation step of organic molecules. [2] Among other things, his research resulted in the famous Sabatier principle, which states that the catalytic effi- ciency is related to an optimal bond strength during chemisorption. The bonds need to be strong enough to actually attach the reactant to the surface while not being so strong that it hinders the formation of the product. [3] The predicted catalytic efficiency of different catalysts are often depicted in a ”volcano plot”, where a schematic version is shown in Figure1.
Au Ag
Pd Pt
Cu Ni
Bond Strength
Catalytic Efficiency
Figure 1: A schematic volcano plot, where different elements have been plotted according the the chemisorption bond strength with a specific molecule and how well the element can catalyse the reaction that the specific molecule will undergo.
However, the bond strength is a quantity that is cumbersome to calculate. It is here useful to instead utilise descriptors, which are quantities that have a high correlation with the property that is to be studied.
In the case of describing the bond strength for a chemisorbed compound on a surface, a descriptor that is widely used is the d-band center of the material. [4] The d-band center in this work is defined to be the median of the d-pDOS. One could even go as far as to again view Figure1, but now picture the x-axis as if it instead represents the d-band center. This might seem like a drastic change, but relation would still to some extent show the same general behaviour. The benefit of considering the d-band center instead of the bond strength is that the d-band center is relatively easy to determine when performing a quantum chemical calculation, and it can even be calculated exclusively for the metal, i.e. it gives an estimate of the reactivity of the metal without having to actually simulate a molecule binding to it. It is also possible to calculate the d-band center for specific atoms in a model, making it possible to get an idea of the reactivity in different regions of the catalyst material.
To experimentally determine and compare the reactivity between different catalytic materials can some- times be difficult. To experientially create a specific nano-particle or surface can be problematic, and the subsequent characterisation can be even harder. Furthermore, when measuring the catalytic efficiently or bond strength during chemisorption the information gained is often an average over the whole sur- face. This is why it is useful to simulate catalytic materials. Here detailed models can be created and studied in order to find trends and develop theories about the nature of the systems. It is when these two approaches are combined in the right way that more efficient research can be achieved. However, there exists a wide variety of different methods that can be utilised in order to simulate a chemical system.
There is a lot of effort in benchmarking and comparing different methods, and there exists methods that are, definitely, more accurate than others, for example the wave function based quantum chemical methods. However, the accuracy often comes at the price of high computational cost, thus leading to long calculation times. Today, in reality, there is no way of calculating a physical-true nano particle with one of these super accurate methods. The reason for this is in large part due to the extreme computer power required to simulate such systems. A widely used computational method is the density functional theory (DFT). It is relatively cheap in computational cost, allowing for the simulation of relatively large systems, and is accurate enough for a range of applications. Still, the computational cost is often too high for calculations of sufficiently large systems and even if a moderate system is to be calculated, one might want the results in seconds or generate them from a regular laptop. In an effort to counter these limitations, a tight binding approximation can be derived from to the original DFT formulation, giving rise to the self-consistent charge density functional based tight binding method (SCC-DFTB).
A more general problem that is encountered when performing quantum chemical calculation is that one generates a lot of data, often making the results hard to interpret. This especially applies when calculating an electronic structure of a system. Looking at the band structure and DOS for Au shown in Figure2, one might get lost in all the curves. Even if relevant information can be read out from this complicated depiction, one might instead choose to find an alternative way of representing the system, while still retaining the important information about the system. Here one again finds use of descriptors.
Conveniently, the d-band center is in close connection to the electronic structure. However, the use of this descriptor, practical as it may be, also raises a few questions. How reliable can a single value be? It might be possible to loose vital information if only a single value is considered. But on the other hand, it might be the case that most other information actually is irrelevant in certain scenarios.
X G W L G
k-points 10
5 0 5 10 15
E (eV)
Density (a.u.)
Figure 2: A representation of the electronic structure of Au. The left side shows the band structure. The right side shows the total density of states (black) and the projected density of states for s (red), p (green) and d (blue).
2 Aim
The aim of this project is to investigate how parameters in the SCC-DFTB method can be tuned to give accurate results for the electronic structure of some transition metals. There will be an attempt to find representative parameters for nickel in a face centred cubic bulk structure, where the reference is cal- culated from DFT. There will also be an effort to understand how these parameters affect the electronic structure and whether or not it is possible to draw physical/chemical conclusions from them. Once work- ing parameters for Nickel has been found, these parameters will be used to simulate larger systems in the form of nano-particles and a surface slab. There will also be an attempt at finding parameters for Copper, Palladium, Silver, Platinum and Gold.
3 Theory
3.1 Molecular orbitals
The forming and breaking of bonds between atoms is something very important in chemistry, one could even state that this is what a chemical process is. But describing what bonds actually are is not trivial.
In an attempt to understand chemical bonding, a theory has been developed called molecular orbital theory (MO-theory). According to MO-Theory, when two atoms are close enough for their individual atomic orbitals to overlap, their orbitals can combine to create molecular orbitals, see Figure3. Two molecular orbitals are formed for each pair of atomic orbitals that combine. Further, the atomic orbitals can combine in two ways. One where the atomic orbitals wave function have the same sign, creating a constructive interference between the atomic orbitals which leads to a bonding molecular orbital. The other possible outcome of the interaction is when the overlapping region of the two atomic orbitals have wave functions of opposite sign, giving rise to destructive interference that results in an anti-bonding molecular orbital. In contrast to the degenerate atomic orbitals, these new bonding and anti-bonding orbitals will have different energies, depicted by their eigenvalues, where the bonding will be lower in energy then the anti-bonding orbital. [5]
Let us consider the valence s-, p and d-atomic-orbitals and the bonds which appears between the same or- bital types for a diatomic molecule. The s-orbitals combine to create a pair of bonding and anti-bonding σ -orbitals. The p-orbitals will combine in pairs of the corresponding atomic orbitals for each atom to create one pair of bonding and anti-bonding σ -orbitals and two pairs of bonding and anti-bonding π- orbitals. The perhaps most important bonding-orbitals for this work is the δ -orbitals which are formed when combining the atomic d-orbitals. The importance of these bonds will be further explored in a later section.
Apart from bonding and anti-bonding orbitals, atomic orbitals can also combine to give rise to non- bonding orbitals. These non-bonding orbitals often refer to the molecular orbitals that has energy eigen- values between the bonding and anti-bonding molecular orbital.
2 3 n
(a)
1
k
Figure 3: (a) An illustration of how an increase in atoms participating in bonding creates more energy eigenstates in the form of molecular orbitals and when combing n orbitals gives rise to n molecular orbitals. The arrows in the figure represents an electron occupying the orbital. The red line shows the highest occupied molecular orbital. (b) The red line depicts a σ - orbital band and corresponding density of states that is formed from combining the s-atomic-orbital of n atoms. The green line represents the σ -orbital band and corresponding density of states that arises from combing the in plane p-orbital of n atoms.
3.2 Band structure
The origins of a band is illustrated in Figure3a. Starting with one atom which is in possession of one electron placed in the s-orbital. By placing an identical atom next to the first one, the two atoms will combine one orbital each to create a bonding orbital and an anti-bonding orbital. The electrons from each atomic orbital will here move to occupy the energetically lowest available molecular orbital. When adding another atom and allowing one orbital from each of the three atoms to combine will give rise to three new orbitals, one bonding, one non-bonding and one anti-bonding, with two electrons occupying the bonding orbital and one electron occupying the non bonding. When n atomic orbitals combine, n molecular orbitals will be created. Here, the energy difference will be small between each eigenstate.
[5]
By introducing lattice points placed in a point of symmetry represented by the atoms, makes it possible to construct a wave function that accounts for the periodicity in the system. This is done by utilizing Blochs theorem, shown in Equation1.
ψk=∑
n
eiknaχn (1)
In the equation, ψk is the wavefunction that is constructred from a linear combination of atomic wave- functions χncentered at lattice point n, separated by the lattice spacing a. k represents a quantum number which determines the phase between different unit cells. In Figure3a the lattice points could be placed at the center of each atom. This placement of the lattice points and by using equation1makes it possible to describe specific phase configurations in terms of ψk. When the overlapping parts of the orbitals all have the same sign then that configuration is in-phase. In the opposite scenario, when all the orbitals overlap with opposing signs, the configuration is out-of-phase. Furthermore, these different phase con- figurations can be viewed as molecular orbitals that are different degrees of bonding, non-bonding and anti-bonding. Where the fully in-phase configuration represents a fully bonding orbital and as the phase shifts the nature of the bond will move towards non-bonding followed by anti-bonding.[6]
When considering Blochs theorem the molecular orbitals instead become crystal orbitals. When dis- cussing crystal orbitals, it is common to refer to the Fermi energy when describing the energy distri- bution of electrons in the crystal orbitals, which is marked by the red line in Figure3a. This is useful because the Fermi energy is the energy of the highest occupied crystal orbital. By instead viewing the electronic structure as crystal orbitals an extension of the model depicted in Figure3a can be achieved, which is shown in Figure3b. In this figure the array of energy eigenstates have been plotted in a region of reciprocal space, known as the first Brillouin zone. Now the quantum number k comes into play, be- cause this number corresponds to a specific coordinate in the first Brillouin zone, which is also referred to as k-space. This leads to each eigenvalue having specific placements in first Brillouin zone. It is when the energy eigenvalues are placed within the context of reciprocal space that one obtains bands. Thus, all the energies that make up the curve results form specific phase configurations, which occupy unique locations in the first Brillouin zone, as illustrated in Figure2.
Studying a band structure more carefully makes it possible to obtain information about the orbital overlap that is involved. The first feature of interest is the dispersion of the band, where a disperse band indicates that a large orbital overlap is present. The second feature of interest is best explained using Figure3b.
In the figure that the red band has a minimum as a starting value and increases in the direction of k while the green band shows the opposite trend. These behaviours are caused by the symmetry of the atomic orbitals. The change in sign of the wave function within p-orbitals causes the pσ orbital to be anti-bonding when the the p-orbitals are aligned after the same sign, which is the case at the center of k-space. The binding pσ orbital is instead found when move away from the center of k-space where every other aligned orbital has a changed sign. The opposite is true for s-orbitals, here the orbital only
consists of one lobe. This means that the orbitals will be in-phase at the center of k-space, and become more out-of-phase as individual orbitals change sign at increasing k. [6]
Moving a step beyond the adding up of atoms in one dimension, toward a three dimensional crystal.
Comparing the schematic band structure in Figure3b to the band structure shown in Figure2, it would not be an understatement to say that the band structure gets a bit more complicated when applied on a real bulk metal. First of all, the labelling of the x-axis has changed from just an increase in k which represented moving away from the center of the first Brillouin zone. This is because there exists many different directions which lead away from the center of the first Brillouin zone. Luckily there also exists high symmetry points that are of interest when sampling the first Brillouin zone for energy eigenvalues.
In order to obtain a proper band structure as the one shown in Figure2, this sampling should be done on lines connecting these points. The most important point of symmetry is the gamma point, depicted by G in all the graphs in this report. The gamma point is extra important since it is this point that represents the center of k-space. Even if the issue of the k-space sampling has been addressed, one might still be a bit uneasy about the sheer complexity of the band structure present in Figure2. Unfortunately this is an effect of moving on to a bulk metal. The bands in these system are not necessary made up from just one type of orbital. The bands in Figure3b are made from s-orbitals combining to form σ -orbitals and p-orbitals combining to form σ -orbitals. But in a metal different kinds of atomic orbital can combine to create crystal orbitals. This of course makes everything much more complicated, but fortunately the curvature can still give information of the nature of the overlap taking place. [6] [7]
3.3 Density of states
The density of states (DOS) is the number of states available for electron occupation at specific energies, and it is derived from the band structure. One can conceptually interpret the DOS as being an alternative way of looking at the band-dispersion where a disperse band gives a low density spread over a large energy range and a flat band gives a high density over a small energy range. The DOS is then compiled of all electronic states. This is shown in Figure3b where the right part shows the DOS. It is also possible to access the contribution of each orbital type in the DOS. This is done though projecting the states on each atomic orbital, giving rise to the conveniently named projected density of states (pDOS). To clarify, in the right part of Figure3b, the area confined by the y-axis and red line would represent the s-pDOS while the area confined by the y-axis and green line would represent the p-pDOS, thus the DOS would be the sum of both areas. The DOS is also of experimental significance, because it is possible to measure the DOS using ultraviolet photoemission spectroscopy (UPS)[8].
By looking at the individual partial density of states makes it possible to get an idea of the orbital-mixing giving rise to individual bands. If for example the s contribution of the density is high at a certain energy, then some of the bands crossing that region should be made out of s-mixed states. This is easy to see in Figure3b and is, though it might be harder to make out, even more useful when trying to understand the electronic structure of Au shown in Figure2.
3.4 The d-band model
The d-band model is motivated by the similarities that the d-states share with atomic orbitals. The bonding that involve atomic d-orbitals in the metal is based on a small overlap, which is shown by the narrow d-pDOS that is created. This causes the d-states to be very local, making it possible to view them as atomic orbitals distributed locally at each atom site in the metal. It is these local orbitals that enables chemisorption of molecules. Forming bonding and anti-bonding molecular orbitals. The width of the d-pDOS, also depend on the placement of the center of the d-pDOS in respect to the Fermi energy, which is illustrated in Figure4. A down shift of the d-band center away from the Fermi energy will cause the
d-pDOS to widen, and in turn making the d-states less local. More delocalised d-states would in turn make chemisorbtion less favorable. [9]
The placement of the center of the d-pDOS in respect to the Fermi energy is called the d-band center, and this single value is central in the d-band model.
Another general behaviour that has been observed is if the d-pDOS ends below the Fermi energy then all the d-states are filled which leads to the chemisorption being more unfavourable. This is due to no orbitals being available for bonding with molecules. If the Fermi energy instead crosses the d-pDOS then there are empty d-states on the surface available, enabling chemisorption of molecules. [10]
Apart from the filling of d-states, the d-band center gives more information. A density of states is made up of bonding, nonbonding and anti-bonding states. As suggested in the previous discussion about the bands, the most bonding states tend to be located at the lowest energies in the DOS, and as the energy increases the less bonding the electronic states will be. This leads to the highest energy states in the DOS being mostly anti-bonding. By placing these states in context with the Fermi energy, one can determine whether the bonding nature of the highest occupied states and also the lowest unoccupied states. If a large portion of the DOS is under the Fermi energy then both the highest occupied and lowest unoccupied states will be anti-bonding making interactions with molecules weaker. In agreement with this, both the highest occupied and lowest unoccupied states will be bonding if the Fermi energy is located at the lower end of the DOS, resulting in more strongly bonded atoms and molecules. [9]
E (eV)
Density
EF
Figure 4: A schematic figure depicting the projected densities, and the effect of shifting the d-band. The Fermi energy (EF) is marked by the green line. The blue curve represents the d-pDOS of the catalyst where the d-band center is shown as the purple line and the red is the s-pDOS of the catalyst. The grey line represents the pDOS of the interacting molecule, where the lower density region is of bonding character and the higher being anti-bonding. The arrows illustrates trend of how the densities change when the d-band center is shifted downward.
The above mentioned properties only take the metal itself into account. But it is also important that these arguments hold when an interacting molecule is introduced. The general relation between the d-pDOS and the orbital of the interacting molecule is illustrated in figure4. Here the grey curve depict the pDOS of the molecule. This pDOS varies depending on the molecule that is interacting. The pDOS will have two high density areas, one bonding and one anti-bonding. When the d-band center is shifted down, the d-pDOS expands downwards becoming more stretched out. This causes the anti-bonding region in the pDOS for the molecule to shift down. Causing the anti-bonding orbitals to be pushed down below the Fermi energy and thus allow for occupation of anti-bonding orbitals. This in turn leads to a decrease in the chemisorption energy. Since these trends are general regardless of the interacting molecule, it is valid to just observe the d-pDOS when evaluating the reactivity of a surface. However, there will still be a difference in the bond strength between different molecules. But the general trend can and should be used as a tool when exploring uncharted interactions. [9]
4 Method
The DFT calculations were performed using the Vienna Ab initio Simulation Package (VASP) [11] man- aged through the Atomic Simulation Environment (ASE) [12]. Similarly, the SCC-DFTB calculations for testing the generated SK-tables were performed using the DFTB+ code [13] and managed through ASE.
A large part of the work has been devoted to create a master script to manage all necessary calculations involved in creating and testing the Slater-Koster (SK) table.
4.1 DFT
DFT is a quantum chemical method that is based on the principle that the total energy can be determined from the electronic density alone [14][15]. Equation2 shows how to derive the energy of the system from the electronic density.
E[n(r)] = Ts+ Eext+ EH+ Exc (2)
In the equation, n(r) is the electronic density. E[n(r)] is the total energy to be obtained, Tsis the kinetic energy, Eectis the external interaction, EHis the Hartree energy and Excis the exchange/correlation term.
The main problem with DFT is to obtain Exc. EDFT would give the exact total energy with the right formulation of Exc, however no such formulation has been found. Instead, there exists a wide variety of approaches (also known as functionals) that are used in order to give a fair representation of Exc. The functional used for the DFT calculations in this report was the generalised gradient approximation Perdew-Burke-Ernzerhof (GGA-PBE) which is a widely used functional. DFT also has shortcomings when predicting unoccupied (virtual) states, which is something that is important to keep in mind. [16]
Two DFT calculations were performed on each system. The fist calculation was done in order to con- verge the electronic structure. Here a Monkhorst pack k-point grid [17] was used for sampling the Brillouin. The second calculation was performed using the previously converged electronic structure and k-sampling on straight lines connecting the special k-points in the k-point grid.
4.2 SCC-DFTB
SCC-DFTB is derived from the total energy expression of DFT by making a Taylor expansion around the ground state electron density. This makes the total energy a functional dependent on (δ n(r)), which is shown in Equation3.
E(δ n(r)) = Eband(δ n(r)) + Ecoul(δ n(r)) + Erep (3) In this equation, Eband(δ n(r)) is the band energy term, Ecoul(δ n(r)) is the additional Coulomb energy caused by the charge fluctuations and Erepis the repulsive energy contribution, which is mainly caused by the repulsive nuclei interactions. [18]
In practice, the different contributions to the total energies are pre-calculated and stored in SK-tables, which speeds up the calculations significantly, making SCC-DFTB faster than DFT. To improve the accuracy of the SCC-DFTB, the construction of these tables allows for some flexibility. Here pseudo atoms are utilised, rather than free atoms, when calculating the relevant integrals. This is done by adding
a confinement potential to the DFT Hamiltonian, to mimic the environment of a bound atom. Typically, a quadratic potential is used, which is shown in equation4.
Vcon f(r) = (r
r0)2 (4)
Here Vcon f(r) is the confined potential which confines the wave function as a function of r. The param- eter to be determined is then r0also known as the compression radius, which then enables the user to influence the confinement. This procedure greatly enhance the accuracy of the method, but at the same time introduce a risk of problems with transferability.
The parameters to be explored in this work are the electronic configuration and the wave function con- finement radius (r0). It is also of interest how these parameters affect the electronic structure of some transition metals. Furthermore, the transferability of these parameters have also been tested by calculat- ing the d-band centres for each atom in Ni metal nano-particles and for atoms in a Ni metal surface with different defects.
4.2.1 Optimisation
Skopt is a software that can optimise the parameters for pairs of atoms. The parameters that are optimised are the compression radii and energy eigenvalues for the s, p, d and f orbitals. These parameters can also be changed manually, but skopt is a tool for automatically optimising these parameters according to some specified rules and criteria. In general, one chooses certain reference-values (in this case generated from the DFT calculation), which is a quantity that can also be obtained through a SCC-DFTB calculation, and the criteria is then that these values should compare in a certain way in order to generate a good fit.
If the value of the quantity that is generated by the SCC-DFTB calculation is dependent on the different compression radii, then there exists a configuration of compression radii or a set of these configurations that will generate the best fit. The optimisation is then an iteration of generating different compression radii configurations so that the best fit is obtained. Skopt can use different optimisation methods.
The method used in the calculations covering this report is particle swarm optimization (PSO) [19]. In this method, a specified number tiles are specified, also known as ”particles”, and how many iterations of optimisation is to be performed. Each particle will over each iteration move towards a local optimal fit. The best global fit will be remembered over all iterations and will be replaced if a better fit is found.
The individual particles will also remember the neighbouring particles fit in order to steer towards the globally best fit.
The rule is then that the compression radii needs to be within specific intervals, which together with the fit makes up the ”multivariable landscape” that the particles are allowed to traverse.
A calculation of the atomic state has to be performed before initiating the skopt procedure. This is done with a software called skgen. It is also here that the electronic occupation is specified. Once the atom has been solved, skopt can be used to generate series of SK-tables.
4.3 Calculating the d-band center
In this report the d-band center is defined as the energy where half of the occupied d-states have lower or equal energy as the energy of the d-band center. This would correspond to an approximate energy of the non-bonding states. The analytical expression is shown in equation5.
2 Z dcent
E0
dPDOS(x)dx = Z EFermi
E0
dPDOS(x)dx (5)
Where E0is the lowest energy in the band-structure, dcentis the d-band center, EFermiis the Fermi energy and dPDOS(x) is the projected density on the d-orbitals as a function of an energy x.
The actual calculations of the d-band center was done in the following way. An initial calculation had to be done by integrating from the lowest energy in the band structure up to the Fermi energy in regard to the partial density belonging to the d-orbitals. By dividing this value by 2 one obtains the density up until the d-center, in accordance equation5. To then determine the d-center one performs an additional integration from the lowest registered energy up to the energy that makes the integral become half of the total density. The upper integration limit is then the d-band center.
5 Results and Discussion
5.1 DFT
The d-band centres from the DFT calculations for Ni, Cu, Pd, Pd, Ag, and Au are given in Table1. The band structure and DOS from the bulk Ni calculation are shown in Figure5. The results in Figure5 and the d-band centres that are tabulated in Table1correlate rather well with the d-pDOS presented by B. Hammer and J.K Nørskov [9]. What makes this even more interesting is that the d-pDOS shown in their article was calculated from metal surface structures. Thinking about the tabulated d-band centres, one might think that Au should have a lower d-center than Ag, which is in contradiction to the values presented. But comparing their d-pDOS one finds that they should both have a complete filling of the d-states. There are actually more factors playing a part in the nobility of metals, for example cohesive energy and matrix coupling [10]. However, this is beyond the scope of this report.
X G W L G
k-points 10.0
7.5 5.0 2.5 0.0 2.5 5.0 7.5
E (eV)
Density (a.u.)
Figure 5: The figure shows the band structure and DOS of Ni calculated using DFT. The coloured lines in the DOS depict the partial projected densities, were red corresponds to s, green to p and blue to d.
Table 1: The energetic position of the d-band center (dcent) gained from DFT calculations. The values correspond to a Fermi energy at 0 eV.
Element dcent(eV)
Ni -1.56
Cu -2.50
Pd -2.36
Ag -4.24
Pt -3.26
Au -3.98
5.2 The influence of the compression radius and electronic configuration on the band A number of calculations on Ni was performed in order to thoroughly examine how the different con- figurations of compression radii and electronic configuration effect the band structure, DOS and pDOS.
In order to achieve an efficient screening, the PSO was constrained so that only a very narrow range of values was allowed from the starting guess.
The calculations that generated the best values for the fit or the best matching d-band-center for the screening calculations on Ni is shown in Table2. The chemical relevance of the d-band center has been stated previously. But one should be careful of reading to much into this one value by itself, and even more caution should be taken if one should choose to model the parameters after the d-band center alone.
Since the d-band center is a value which is supposed to model the center of the d-pDOS, it is possible to obtain the same value from different density of states. Thus obtaining the same d-band center with completely different electronic structures and loosing a deeper chemical relevancy in the results. This problem is apparent when comparing the fit and d-band center for the SCC-DFTB calculations in Table2.
Here a bad fit can still generate a good approximation of the d-band center. It is also the case that a good fit can give a bad approximation of the d-center, ideally a perfect fit would also lead to a perfect match for the d-band center. The fact that no such match has been found indicates that no perfect parameters
have been found, if such a match is possible. However as previously discussed, one can argue that even if the underlying electronic structure might be wrong, it might still be used in a relativistic way, where only the shift of the d-band center is considered.
Table 2: The calculations that generated the best fit or best approximated d-center, dcent. The electronic configurations, ec, used in each calculation is also stated. The naming of the calculations also show the compression radii used and is on the form ”element r0s-r0p-rd0”.
calculation Fit (arbitrary units) dcent(eV) ec
Ni 3-3-6 1.18 -1.58 (Ar)3d84s2
Ni 3-3-4 1.19 -1.44 (Ar)3d84s2
Ni 3-3.7-3 1.07 -1.64 (Ar)3d94s1
Ni 3.5-3.5-3 0.979 -1.70 (Ar)3d94s1
Ni 3.3-3.7-3 1.01 -1.68 (Ar)3d94s1
Ni 2-2-3 4.80 -1.47 (Ar)3d94s1
Ni 3.4-3.5-3 0.951 -1.70 (Ar)3d94s1
Ni 4-4-3 1.08 -1.75 (Ar)3d94s1
Ni 2.5-2.5-8 1.61 -1.57 (Ar)3d84s2
Ni 7-7-3 3.42 -1.58 (Ar)3d84s2
Ni 2-2-8 2.27 -1.51 (Ar)3d84s2
The best fit, Fit= 1.01, was obtained with the compression radius configuration r0s = 3.3, r0p= 3.7 and r0d= 3 together with the electronic configuration (Ar)3d94s1, which generated a d-band center of −1.68.
A comparison of the DOS from this calculation (red line) with the DOS from the DFT calculation (black line) is presented in Figure6a. The band structure generated from this calculation is shown in Figure6b.
Here the colours indicate how well each increment of the band, match with the band generated from the DFT calculation. A blue colour means that the value in the presented band is 1 eV lower than the DFT representation, a red colour means that the energy is 1 eV too high. The closer to the correct value the increment of the band has the darker the colour, where a black colour indicates a small error.
10 8 6 4 2 0 2
E (eV) 0.0
2.5 5.0 7.5 10.0 12.5 15.0
Density (a.u.)
(a)
X G W L G
k-points 8
6 4 2 0 2 4 6
E (eV)
(b)
Figure 6: The accuracy of the electronic structure has been presented in two ways. (a) The DOS generated from the SCC- DFTB calculation (red) and the DOS from the DFT calculation (Black). The vertical lines shows the d-band center for each calculation, red for SCC-DFTB and dotted black for DFT. (b) The SCC-DFTB generated band structure, were the error in comparison to the DFT generated band structure is illustrated by the colour. Black indicates a match between the two band structures. Red indicates that the energies in the k-point is to high in comparison to the reference and blue indicates that the energy is too low.
There was one set of parameters that generated what looked to be an even better matching band structure, see Figure7. This was done by using configuration r0s = 3.3, r0p= 3.8 and r0d= 3.3, which resulted in a Fit=0.84853 and dcent= −1.7624.
10 8 6 4 2 0 2 E (eV)
0.0 2.5 5.0 7.5 10.0 12.5 15.0
Density (a.u.)
(a)
X G W L G
k-points 8
6 4 2 0 2 4 6
E (eV)
(b)
Figure 7: The accuracy of the electronic structure has been presented in two ways. (a) The DOS generated from the SCC- DFTB calculation (red) and the DOS from the DFT calculation (Black). The vertical lines shows the d-band center for each calculation, red for SCC-DFTB and dotted black for DFT. (b) The SCC-DFTB generated band structure, were the error in comparison to the SCC-DFTB generated band structure is illustrated by the colour. Black indicates a match between the two band structures. Red indicates that the energies in the k-point is to high in comparison to the reference and blue indicates that the energy is too low.
Comparing Figure6and7, the latter does look to better represent the electronic structure on average.
However, the fit was better for the calculation that gave rise to the electronic structure depicted in figure 6. This is a consequence of the criteria for the fitting, and how well these criteria translate to what is visually interpreted as a good match. The average error might be lower while having outlier values, which are not heavily weighted, thus still generating a good fit. It was decided an even distribution of the errors is preferred since this was thought to give a better representation of the overall electronic structure. However, if Fit=1 exactly, then every value in all the bands must match perfectly and as long as this is not the case the fit should be evaluated carefully. It is thus important to decide whether a band structure with a good fit really is the best representation of the reference band structure. This is more important since the fit is only dependent on direct values while the curvature of the bands (derivatives) are of chemical importance. By further comparing the DOS for these calculations one can see that the best fitted calculation also gives a better d-band center at the cost of not approximating the lower end of the DOS very well. One might argue that the lower states will not participate much in interactions and are not that chemically important. But if the goal is to create a physically representative electronic structure an easy feature to evaluate would be the energy span of the DOS, hence making the lowest and highest bands in the DOS important. Further, while the d-band center is of chemical importance, it is the trend in the shifting of the d-band center between the bulk and different surfaces or nano-particles that are important. And the hope in the parametrisation is that the individual orbitals can fundamentally be represented using the same compression radii regardless of environment. That is, provided the representative compression radii for the pairwise interaction, a correct electronic structure will emerge for the specific environments.
The resulting band structure of Ni ((Ar)3d94s1) for different r0is plotted together with the error, which is represented by the colour, is shown Figure8a-d. Where Figure8b shows the band structure when rs0 is changed,8c shows the band structure when r0pis changed and8d shows the band structure for a varied r0d.
X G W L G k-points
8 6 4 2 0 2 4 6
E (eV)
(a)
X G W L G
k-points 8
6 4 2 0 2 4 6
E (eV)
(b)
X G W L G
k-points 8
6 4 2 0 2 4 6
E (eV)
(c)
X G W L G
k-points 8
6 4 2 0 2 4 6
E (eV)
(d)
X G W L G
k-points 8
6 4 2 0 2 4 6
E (eV)
(e)
Figure 8: The different band structures in the figure illustrates how varying each parameter alters the band structure. (a) Is the starting configuration, using the electronic structure (Ar)3d94s1and a compression radii configuration of r0s≈3, r0p≈3 and r0d≈3. (b) r0s was changed to ≈ 3.3 . (c) r0pwas changed to ≈ 3.3. (d) r0dwas changed to ≈ 3.3. (e) The electronic structure was changed to (Ar)3d84s2.
The band depicted in Figure8e uses the same r0configuration as in Figure8a, but was created with the electronic configuration (Ar)3d84s2. One can clearly see that the energies tend to shift upward when the configuration (Ar)3d84s2is applied.
Comparing the band structures in Figure8, it appears that the different compression radii r0s, r0pand r0dare not independent of each other. Similar features in the band structure are altered when these parameters are varied individually. However, there are some subtle differences in how the different parameters alter the band structure. rd0seems to have a large influence on the framework for the band structure, where the energy span between the individual bands are set. r0s has much influence on the dispersion of the lower bands and r0phas much influence on the dispersion of the bands around the Fermi energy. Another way of describing these tendencies is to state that r0d has most weight in describing the energetic positions.
While r0s, at low energies, and r0p, at higher energies, has most weight on describing the derivatives of
these energetic positions. But as stated above, the parameters are not completely independent. Although, the different compression radii can change the energetic positions and the dispersions, they do so at different magnitudes. When r0dis increased the bands get shifted downwards to lower energies, while at the same time making the lowest bands more disperse and the bands above less disperse. An increase of r0s shifts the lowest bands upwards making them less disperse and also shifts some of the higher bands downwards. When increasing r0psome bands around the Fermi energy gets shifted up while also making the bands in the middle more disperse.
The above stated behaviours might be a bit puzzling when viewed in a chemical context. We address r0s in the context of the band structure and interpret it from a chemist’s point of view. For Ni with the atomic configuration (Ar)3d94s1 it is reasonable to assume that the lowest energy interaction is a σ - bond. Indeed the lowest band in the DFT band structures for Ni is a disperse band with a minimum in the gamma point which indicates a bonding σ . If this band represents the phase of the s orbitals, then one might expect the band to be even more disperse when the overlap is increased (Figure8a and b). But this is not what is observed, instead the band gets less disperse. The reason for this might be that too much overlap between orbitals creates an interaction that is less like a σ -bond and that is less favourable, thus making the bonding interaction less bonding. According to the electronic configuration, there should not be any occupation of p orbitals. Nonetheless, r0pstill have an effect on the electronic structure. There is also p-state contribution in the DOS for the DFT calculation (Figure5). This might be due to a relaxation that by p-mixed interactions that take place when the atoms are placed within the crystal structure. From the changes between Figure8a and c one can also deduce that these p-mixed interactions produce the disperse bands stretching from ≈ 2 eV to the Fermi energy and also the bands above the Fermi energy.
Furthermore, the bands under the Fermi energy are anti-bonding in gamma and the ones above bonding, which further indicates that π-bonding is involved. That r0d influences the energetic positions of most of the bands is a bit counter-intuitive, if one studies the projected density of states in Figure5. The d- pDOS covers large parts of the energy span of most of the occupied bands, which might indicate that the covered bands have a large contribution from d-orbitals. However according to the d-band model, the d-states should be very local and behave as atomic orbitals. Viewing the down shift of the bands when rd0 is increased (figure8a and d) in accordance with the d-band model. One might speculate that the down shift in energy is due to de-localisation of the d-orbitals, allowing for further mixing of states which gives rise to a relaxation of the electronic structure. If this is the case then the r0dcould be considered as the most important parameter to get right. Because it is the reactivity of these localised d-states that are important when studying heterogeneous catalysis.
Comparing the band structures with different electronic configurations, Figure 8a and e. The elec- tronic ground state configuration for a free Ni atom is (Ar)3d84s2, while in a solid the ground state is (Ar)3d94s1. As observed the energies should increase in the overall band structure for (Ar)3d84s2 since this is not the favoured configuration of metallic Ni. The fact that the energies are lower for the (Ar)3d94s1 configuration indicates that this is the electronic configuration that is present in metallic nickel.
Large parts of the discussion up until this point is based on one, still, big assumption. That is that the parameters are physically relevant. The tight-binding approximation is just that, an approximation.
To make things worse, it is an approximation of an already approximate method, DFT. Since DFT is a widely used method for describing chemical systems one could perhaps feel comfortable with asso- ciating the results from this method with a physical meaning. But it is possible that the tight-binding approximation might to some extent make the results diverge away from this real context. It is also possible that the orbital concept is only relevant if the right parameters are present, and if the wrong parameters are applied the resulting orbitals might be non-physical. One thing that would help the cred- itability of SCC-DFTB calculations is if one can exactly replicate the results from DFT calculations with a SCC-DFTB calculation. Even if this might be nearly impossible to do, it is still not enough. There
are some more criteria regarding the parameters that also needs to be fulfilled. Firstly only one set of parameters and electronic configuration should be able to give a perfect replication of the DFT results. If more exists then the physical meaning behind the r0configuration would become distorted. If the system has orbitals which can be confined, there should also only exist one set of r0configuration which would lead to a physically real electronic structure. The counter argument to this is that there exist several r0 configuration that leads to this electronic structure, where one is a real solution and where the rest gives the correct answer out of mere coincidence. However, this would only make the method more unreliable when trying to associate chemical properties with the parameters. It would be impossible from the cal- culation alone to decide which r0configuration is correct. Thus, there would always be a risk of making chemical conclusions from non-physical data. Luckily, it seems unlikely to replicate the DFT electronic structure using the (Ar)3d84s2configuration. From the calculations done it also seems that there exists convergence point toward generating a perfect fit against the DFT band structure. But regardless of the physical relevancy of the SCC-DFTB parameters one can still be pragmatic and only view it as a tool that can be used to make calculations easier, however one should then be careful of reading out to much from the results if there doesn’t exist a proper reference.
5.3 Finding parameters using PSO
A performance test of the PSO through Skopt was done in order to determine how appropriate this method was for finding the right parameters. This was done by selecting a small optimisation window where a good fit existed. The PSO was allowed to place 25 particles in the range rs0= 3.3 ± 0.5, r0p= 3.8 ± 0.5 and r0d= 3.3 ± 0.5, then letting the PSO look for the best fit over 10 iterations. This generated a fit of 0.82955. However, this fit indicates that PSO might not be the best approach to find the compression radii. Even when provided with a small operation window where a good fit is located the optimisation generated a worse fit then the manual optimisation. This was the case when using 25 particles over 10 generations, perhaps more particles and iterations are needed in order to find a sufficient fit. In order to make this feasible, one would have to use several computer cores when performing the calculations. It is also important to note that the PSO was performed while given well motivated starting guesses to help find a good fit. If this is not known when setting up the PSO, then it would require a bigger optimisation window resulting in an even harder task for the computer. A better approach might be to initially find a good fit by manually vary the r0configuration and then use the PSO to do the last adjustments in order to generate a next to perfect fit.
5.4 Extending the parametrisation on other elements
Once a sufficient fit for Ni was found, an effort was made to find a way of translating the parameters used to the selection of parameters for Cu, Pd, Ag, Pt and Au. This was done by applying scaled versions of the working parameters for Ni (r0s = 3.3, r0p= 3.8 and rd0= 3.3) on the new elements. The scaling was performed by writing an algorithm that created permutations of scaling factors ranging from 1.14 to 1.38. The best matching electronic structures are presented in Figure9. The parameters that was used in order to generate these DOS are presented for each element in Table3. The table also shows the obtained d-band center and the difference in d-band centre between the SCC-DFTB and DFT calculation for each element. Studying the resulting d-band center one can see that there seems to be a bias present. While the overall DOS matches rather well between the DFT and SCC-DFTB generated DOS, the d-band center is always just a bit lower, which is an effect of how the DOS matched between the two different methods.
It was rather difficult to get all the individual peaks to line up exactly. In all the DOS depicted, the DOS is overestimated at low energies, while the peaks in the middle are shifted down in energy and the high energy peaks are lower in density than the DFT generated DOS. This, of course, would generate a lower d-band center. All this indicates an over-relaxation of the electronic structure, but the question is if this is
something inherent from the tight binding approximation or if this is due to the wrong parameters being applied. If it is something caused by applying the tight binding approximation, then one would have too keep this in mind in further parametrisation, where some eigenvalues could be manually shifted instead.
10 8 6 4 2 0 2
E (eV) 0.0
2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0
Density (a.u.)
(a)
8 6 4 2 0 2
E (eV) 0
2 4 6 8 10 12
Density (a.u.)
(b)
8 6 4 2 0 2
E (eV) 0.0
2.5 5.0 7.5 10.0 12.5 15.0 17.5
Density (a.u.)
(c)
12 10 8 6 4 2 0 2
E (eV) 0
2 4 6 8 10
Density (a.u.)
(d)
10 8 6 4 2 0 2
E (eV) 0
2 4 6 8 10 12
Density (a.u.)
(e)
Figure 9: The density of states for Cu (a), Pd (b), Ag (c), Pt (d) and Au (e). The black line represents the DOS from a DFT calculation while the red line is from the dftb calculation. The red vertical line shows the d-band center from the dftb calculation and the dashed black vertical line represents the d-band center from DFT.
Table 3: The compression radii that generated the best matching electronic structure for the calculated elements. The obtained d-band center dcent and the difference between the SCC-DFTB and DFT d-band center ∆dcentre f is also shown for each calculated element. A negative ∆dcentre f indicates that the SCC-DFTB generated d-band center is lower then the DFT d-band center.
Element r0s r0p rd0 dcent ∆dcentre f
Cu 3.80 4.18 3.63 -2.82 -0.312
Pd 3.96 4.65 3.96 -2.62 -0.261
Ag 4.16 4.65 4.04 -4.45 -0.208
Pt 3.96 4.56 4.16 -3.71 -0.451
Au 4.16 4.79 4.16 -4.36 -0.381
