Linköping University | Department of Physics, Chemistry and Biology Bachelors thesis, 16 hp | Educational Program: Chemical analysis engineering Spring term 2020 | LITH-IFM-G-EX—20/3847--SE
A model for heterogenic catalytic
conversion of carbon dioxide to
methanol
Elin Johannesson
Examinator: Lars Ojamäe Tutor: Karl Rönnby
Division, department
Department of Physics, Chemistry and Biology Linköping University
ISBN
ISRN: LITH-IFM-G-EX--20/3847--SE
_________________________________________________________________
Serietitel och serienummer ISSN
Title of series, numbering ______________________________
Språk/Language Svenska/Swedish Engelska/English ________________ Rapporttyp Report category Licentiatavhandling Examensarbete C-uppsats D-uppsats Övrig rapport _____________ Title
A model for heterogenic catalytic conversion of carbon dioxide to methanol
Author
Elin Johannesson
Keywords
Carbon dioxide, methanol, heterogeneous catalysis, zinc oxide, computational chemistry, density functional theory.
Abstract
Since our society became industrialised, the levels of carbon dioxide in our atmosphere have been steadily rising, to the point where it in early 2020 at is 413 ppm. The high concentration is causing several troubling effects worldwide because of the increase in mean temperature that it creates, which causes longer draughts, more severe floods, and rising seawater levels to name a few. There are a few measures that can be taken to reduce carbon dioxide in the atmosphere, among which there are a number of methods that currently are being researched and/or used. The prospect of capturing carbon dioxide and using it as a carbon building block to make methanol is one solution that is particularly interesting, since it in theory could provide a fuel for combustion engines that is net neutral regarding carbon emission. Methanol can be synthesised from carbon dioxide using a heterogeneous catalyst consisting of copper, Cu, and zinc oxide, ZnO. This research is focused on one of the components of the catalyst, the metal oxide ZnO in the form of crystallites or nanoparticles (ZnO)n.
Quantum chemistry is a branch of computational chemistry which is centered on solving the Schrödinger equation for molecular systems. Density functional theory, DFT, is an approach to quantum theory which in this study was used to calculate the geometry and energy of the particles. The supercomputer Tetralith in the National Supercomputer Centre, NSC, was used to carry out the calculations. The DFT calculations utilized the functional B3LYP and the basis set 6-31G (d,p). One of the largest particle sizes studied, (ZnO)20, with a structure
that has a large, flat surface, was found to be the most energetically favourable. According to studies, the presence of an oxygen vacancy on the surface of ZnO reduces the amount of activation energy required for CO2 to bond to the particle, which increases the chance of
forming CO and thus continuing the process of forming methanol. Two structures of (ZnO)20 were investigated in this regard, where
oxygen atoms were removed at different locations, creating four versions of Zn20O19 in total. This proved yet again that the version with a
large, flat surface yields the lesser amount of energy when an O atom is removed from the centre of its surface. The adsorption of CO2 to
the ZnO clusters was studied by calculating the energy of adsorption, and this showed that it was the second version of (ZnO)20, without
an O vacancy, that yielded the least amount of energy, thus being the most favourable species to engage in physisorption with CO2.
Lastly, the activation energy was investigated, and a diagram of the reaction process of CO2 adsorbing to Zn20O19 forming (ZnO)20 and CO
is presented in this paper, which shows that the required activation energy is 127 kJ/mol.
Date 2020-06-10
Abstract
Since our society became industrialised, the levels of carbon dioxide in our atmosphere have been steadily rising, to the point where it in early 2020 at is 413 ppm. The high concentration is causing several troubling effects worldwide because of the increase in mean temperature that it creates, which causes longer draughts, more severe floods, and rising seawater levels to name a few. There are a few measures that can be taken to reduce carbon dioxide in the atmosphere, among which there are a number of methods that currently are being researched and/or used. The prospect of capturing carbon dioxide and using it as a carbon building block to make methanol is one solution that is particularly interesting, since it in theory could provide a fuel for combustion engines that is net neutral regarding carbon emission. Methanol can be synthesised from carbon dioxide using a heterogeneous catalyst consisting of copper, Cu, and zinc oxide, ZnO. This research is focused on one of the components of the catalyst, the metal oxide ZnO in the form of crystallites or nanoparticles (ZnO)n. Quantum chemistry is a branch of computational chemistry which is centered on solving the Schrödinger equation for molecular systems. Density functional theory, DFT, is an approach to quantum theory which in this study was used to calculate the geometry and energy of the particles. The supercomputer Tetralith in the National Supercomputer Centre, NSC, was used to carry out the calculations. The DFT calculations utilized the functional B3LYP and the basis set 6-31G (d,p). One of the largest particle sizes studied, (ZnO)20, with a structure that has a large, flat surface, was found to be the most energetically favourable. According to studies, the presence of an oxygen vacancy on the surface of ZnO reduces the amount of activation energy required for CO2 to bond to the particle, which increases the chance of forming CO and thus continuing the process of forming methanol. Two structures of (ZnO)20 were investigated in this regard, where oxygen atoms were removed at different locations, creating four versions of Zn20O19 in total. This proved yet again that the version with a large, flat surface yields the lesser amount of energy when an O atom is removed from the centre of its surface. The adsorption of CO2 to the ZnO clusters was studied by calculating the energy of adsorption, and this showed that it was the second version of (ZnO)20, without an O vacancy, that yielded the least amount of energy, thus being the most favourable species to engage in
physisorption with CO2. Lastly, the activation energy was investigated, and a diagram of the reaction process of CO2 adsorbing to Zn20O19 forming (ZnO)20 and CO is presented in this paper, which shows that the required activation energy is 127 kJ/mol.
Abbreviations CO2 – Carbon dioxide CO – Carbon monoxide CH3OH – Methanol ZnO – Zinc oxide HCO3- - Bicarbonate AO – Atomic orbital MO – Molecular orbital
HOMO – Highest occupied molecular orbital LUMO – Lowest unoccupied molecular orbital HF – Hartree-Fock
SCF – Self-consistent field
HF-SCF – Hartree-Fock self-consistent field LCAO – Linear combinations of atomic orbitals CSF – Configuration state function
Table of content
1. Introduction ... 1
1.1 Background ... 1
1.2 Aim... 2
2. Theory ... 2
2.1 Methanol synthesis from carbon dioxide ... 2
2.2 Heterogeneous catalysis ... 2
2.2.1 The Langmuir-Hinshelwood mechanism ... 3
2.2.2 The Eley-Rideal mechanism ... 3
2.3 Methanol synthesis by heterogeneous catalysis... 4
3. Computational chemistry ... 4
3.1 The core of computational chemistry ... 4
3.1.1 The Schrödinger equation ... 4
3.1.2 Molecular orbitals- Linear combinations of atomic orbitals ... 7
3.1.3 Hartree-Fock theory ... 7
3.2 Computational methods ... 10
3.2.1 First approach – ab initio methods... 10
3.2.2 Second approach - density functional theory ... 11
4. Method ... 12
4.1 Computer and programs ... 12
4.2 ZnO-particles: size and structure ... 12
4.3 ZnO-particles: oxygen vacancy ... 13
4.4 CO2 bonding to ZnO ... 13
4.5 Activation energy ... 14
5. Results ... 14
5.1 Size and structure of ZnO ... 14
5.2 Oxygen vacancy ... 16
5.3 CO2 adsorption ... 17
5.4 Conversion of CO2 to CO ... 18
6. Conclusion and discussion ... 20
7. Acknowledgements ... 22
References ... 23
1
1. Introduction
1.1 Background
Carbon dioxide, CO2, is a natural component in the atmosphere of the earth. The molecule is a crucial element in the photosynthesis, and is emitted from animals through respiration, volcanic eruptions, and other natural sources. CO2 is a so-called greenhouse gas, meaning that it has the ability to reflect heat-radiation back to the earth, instead of letting it pass out through the atmosphere. The
“greenhouse-effect” is contributing to making the planet habitable with temperatures that can support life, but the high abundance of carbon dioxide in the atmosphere today is having several troubling consequences [1].
Prior to the industrialisation of our society, the concentration of carbon dioxide in the atmosphere was about 280 ppm [2]. With industrialisation came the burning of fossil fuels, coal, oil, and natural gas, all with high carbon content. In modern society we can also see high CO2 emissions from large industries, deforestation, and the production of cement, to name a few examples beside the non-stop burning of fossil fuels [1]. All these emissions contribute to the steadily rising concentration of carbon dioxide in our atmosphere, which in early 2020 is at a level of 413 ppm [3].
A significant increase of the levels of greenhouse gases in the atmosphere, such as we have seen in the past centuries, can cause the mean temperature of the globe to rise. Rising the mean
temperature, even by only 2°C, will lead to rising seawater levels, an increased risk of flooding, an increase in heatwaves, less fresh water, reduced crop yields, more frequent droughts, and so on. In addition, these consequences will be more prominent in areas around the equator, such as Africa, Asia, and South America, making this not only a problem regarding social aspects, but ethical aspects as well [4]. The ways in which nature itself can control the amount of CO2 in the atmosphere are no longer enough, making human intervention necessary.
There are three main tactics when it comes to reducing the CO2 concentration and emissions. One is to capture and store CO2 to lower the concentration. Another way is to develop and improve alternative sources for fuel such as hydrogen for example, to decrease the emissions. A third option is to capture and use the CO2 as a carbon building block, for an example in the synthesis of more useful compounds such as methanol, dimethyl ether, or hydrocarbons. If methanol is synthesised using captured CO2 and hydrogen that has been produced using a renewable source of energy, it can be used as an alternative fuel which then would be net neutral in regard to CO2 emissions [2]. Commercial methanol is generally produced using the method of synthesis gas (“syngas”), where a gas mixture of CO, CO2 and H2 is hydrogenated using heterogeneous catalysis where the catalysts consist of copper, Cu, and zinc oxide, ZnO, at 5-10 MPa and 200-300 degrees Celsius. The synthesis gas-method and its reaction has been well studied and optimised to fit the needs and demands of the industry and its economic factors. When only using CO2 and H2 to synthesise methanol however, there are different requirements and the mechanism of the reaction is different. Thus, the synthesis of methanol using CO2 demands further research and optimisation to make it commercially,
2
1.2 Aim
The synthesis of methanol from CO2 using heterogeneous catalysis is a reaction which has been well studied in recent years, but some of the fundamentals of the reaction are still relatively unknown. The aim of this study is to research the most energy-sufficient way to synthesise methanol from carbon dioxide and hydrogen using heterogeneous catalysis, with a focus on the catalytic component zinc oxide, ZnO. The study will therefore focus on researching the more energy-sufficient size and structure of ZnO-particles, whether an oxygen vacancy in the cluster is energetically favourable when adsorbing CO2, and to calculate the activation energy required for CO2 to bond with a ZnO-cluster.
2. Theory
2.1 Methanol synthesis from carbon dioxide
From a stoichiometric standpoint, one way to synthesise methanol from carbon dioxide is shown below:
𝐶𝑂2 (𝑔)+ 3 𝐻2 (𝑔)→ 𝐶𝐻3𝑂𝐻(𝑔)+ 𝐻2𝑂(𝑔)
This reaction has, at 298 K, a net enthalpy value of ΔH0= -49,5 kJ/mol, meaning that the reaction is exothermic [2]. The entropy value of the reaction however is ΔSm0= -177 J/K*mol, indicating that even though the reaction is exothermic, the reaction is not spontaneous at 298 K [6]. CO2 has a bond dissociation energy of 1072 kJ/mol, meaning that the required energy input is large [7]. A
consequence of the high bond dissociation energy is that the reaction has a high activation energy, making it a slow reaction. To make this reaction occur one could, from a thermodynamic standpoint, increase the pressure or decrease the temperature in the reaction chamber. This could benefit the synthesis of methanol but comes with an increased risk of forming undesired compounds such as CO, higher alcohols, or hydrocarbons [2].
There are multiple ways to change the rate of a chemical reaction. From studying the Arrhenius equation, equation (1), it is evident that increasing the temperature T or lowering the activation energy Ea will increase the rate constant k of any given reaction [8].
ln 𝑘 = ln 𝐴 − 𝐸𝑎
𝑅𝑇 (1)
A more efficient way to synthesise methanol is through the hydrogenation of CO2 using
homogeneous or heterogeneous catalysts. The homogeneous catalysis, in which the catalyst and the reactant are in a similar phase, is functional in this aspect but the heterogeneous catalysis is
preferred because of its stability, regeneration and lower costs regarding large scale production [2]. By definition, a catalyst is a substance that can enable a reaction to occur at a higher rate or under different conditions that is usually required. The catalyst is a part of the reaction, but it will not be consumed or changed in any way, meaning that once the reaction has stopped there will have been no net change for the catalyst [8]. By introducing a catalyst in a chemical reaction, the reaction is provided a different path and by that the activation energy Ea that is required for the reaction to take place is lowered. As a result, the rate of the reaction will increase [8].
2.2 Heterogeneous catalysis
A heterogeneous catalysis is when the catalyst used in a reaction has a different phase from the reactants, solvents and other substances that are parts of the reaction. In this study, it is the reaction
3
between gaseous molecules using metallic or metal oxide surfaces as a catalyst that is relevant and will be focused on.
The mechanism by which this type of catalyst interact with the reactants in a reaction is by
adsorption. The reactant, at least one of them but sometimes both, is adsorbed to the surface of the metal, where it usually is modified to an active form which in turn can react with the other reactant thus forming the product and will subsequently be desorbed from the surface [9].
There are two different ways in which atoms and molecules can attach themselves to for example a metallic surface. Physisorption, short for “physical adsorption”, is when the atom or molecule can participate in van der Waals interactions, e.g. dipolar interactions or London dispersion, with the metal. On the positive side, these types of interactions have a long range which sometimes is necessary for a reaction to occur. Unfortunately, the van der Waal interactions are weak.
Physisorption is not desired in the reaction where carbon dioxide and hydrogen forms methanol, because the small enthalpy changes of the van der Waal bonds are not sufficient to break the intramolecular bonds required [9].
In chemisorption, short for “chemical adsorption”, on the other hand, the atom or molecule creates a chemical bond with the surface. Covalent bonds are the most common chemical bonds involved, and this contributes to the fact that the distance between the atom/molecule and the surface is usually much shorter than in physisorption. These types of bonds are also much higher in energy than the van der Waals interactions of physisorption, the typical enthalpy values for chemisorption are about 10 times higher than the highest observed enthalpies for physisorption (20 kJ/mol physisorption, ~-200 kJ/mol chemisorption). The high enthalpy values of these bonds in chemisorption is the main reason as to why molecules can attach to the surface of the metal and then break up into fragments, which in turn allows them to react with other reactants and eventually forming a product [9]. There are two main mechanisms of reactions that are surface-catalysed, these are called the Langmuir-Hinshelwood and the Eley-Rideal.
2.2.1 The Langmuir-Hinshelwood mechanism
The Langmuir-Hinshelwood mechanism is when two reactants, A and B, has adsorbed to the surface, S, of the catalyst. The reactants initially form the adsorbed substances AS and BS, and the reaction to create the product, P, occurs when AS and BS encounters [9].
𝐴 + 𝑆 → 𝐴𝑆 𝐵 + 𝑆 → 𝐵𝑆 𝐴𝑆 + 𝐵𝑆 → 𝑃
2.2.2 The Eley-Rideal mechanism
The Eley-Rideal mechanism takes place when only one of the reactants A and B has adsorbed to the surface, S. The remaining non-adsorbed reactant then reacts with the adsorbed reactant and forms the product, P [9].
𝐴 + 𝑆 → 𝐴𝑆 𝐴𝑆 + 𝐵 → 𝑃
4
2.3 Methanol synthesis by heterogeneous catalysis
As previously mentioned in 1.1, the synthesis gas reaction of CO/CO2/H2 to methanol using heterogeneous catalysis is well studied. When only using CO2 and H2 to synthesise methanol however, the requirements for the reaction and the reaction paths are different. The research and optimisation of this method is ongoing and is generally focused on active sites, particle size/structure and the interface between Cu and ZnO [5].
Regarding the optimisation of catalysts, there have been multiple metals that has been researched through the years. Today, Cu is commonly regarded as the most suitable metal to use as an active catalyst component, but it is normally used with added modifiers such as Zn, Al, Ga, and Si to name a few. The modifiers are usually in the form of oxides and they are used to help form and stabilize the active phase of the catalyst and will also help modifying the interactions that occur between Cu and the reactant. In zinc oxide, ZnO, vacancies have been observed in its structural lattice from oxygen, leaving an electron pair which then can interact with reactants. ZnO can also stabilize Cu and enhance the dispersion [2]. Dispersion is a term that describes the ratio between the number of surface atoms and the total number of atoms. In a small molecule, for an example in the size of three or four atoms, all atoms are surface atoms. When studying larger molecules, some atoms are bound to be surrounded by adjacent atoms and will therefore no longer be on the surface of the particle. The dispersion can also depend on how the atoms are packed, and of the shape of the particle [10].
𝒟 =𝑛𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑠𝑢𝑟𝑓𝑎𝑐𝑒 𝑎𝑡𝑜𝑚𝑠 𝑡𝑜𝑡𝑎𝑙 𝑛𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑎𝑡𝑜𝑚𝑠
There are a few uncertainties regarding the mechanism of the reaction in which methanol is synthesised from CO2 using heterogeneous catalysis. The part of the reaction that is currently accepted however, is that the hydrogen is adsorbed and dissociated on the Cu surface, while CO2 is chemisorbed to ZnO and forms bicarbonate, HCO3-. Atomic hydrogen is then transferred from Cu to ZnO, where it hydrogenates HCO3- to eventually form methanol, CH3OH. To fully understand the reaction that occurs and through that being able to optimise this reaction, extensive research is continuously carried out with a focus on the interface between Cu and ZnO, and also identifying active sites on both species [7].
3. Computational chemistry
3.1 The core of computational chemistry
Increasingly intricate and complex questions regarding for an example molecular geometry, energy, and reactivity, requires a more complicated and comprehensive approach. A collection of techniques that makes it possible to study chemical problems on a computer, called computational chemistry, makes this possible. Computational chemistry is common in chemical fields such as nanotechnology, materials science, and pharmaceuticals to name a few.
There are different groups of techniques within the field of computational chemistry. The perhaps most relevant of them is quantum chemistry, which can be divided into ab initio calculations,
semiempirical calculations and density functional theory calculations. At the centre of these different methods of quantum chemistry lies the Schrödinger equation [11].
3.1.1 The Schrödinger equation
While classical physic treats electrons as particles, several experiments in the early 1900´s changed the way science perceived matter and has later been considered the starting point of modern-day physics. In 1925, the Davisson-Germer experiment was carried out, in which the diffraction of
5
electrons by a crystal was observed, indicating that particles have wave-like characteristics. The conclusion was that particles has properties which in classical physics are attributed to waves, and that electromagnetic radiation has some properties classically viewed as attributes to particles, thus the wave-particle duality was formed [12].
Quantum mechanics was created to take place where the classic physics fell short, namely in
microscopic systems such as molecules or atoms with small masses. Erwin Schrödinger presented the Schrödinger equation in 1926, proposing that it could solve the wavefunction of any atom or
molecule. The time-independent Schrödinger equation in one dimension, where the system remains unchanged over time, is presented in equation (2):
− ħ 2 2𝑚 𝑑2𝜓 𝑑𝑥2+ 𝑉(𝑥)𝜓 = 𝐸𝜓 (2) 𝐸 is the total sum of the kinetic and potential energy of the system, and 𝜓 represents the
wavefunction. Since Etot= Ek+Ep, it can be deduced that the first term in the equation is the kinetic energy, and 𝑉(𝑥) is the potential energy of the system at a point “x”. The dynamics of a system, such as momentum and location, is described by the wavefunction. Regarding location, the Born
interpretation states that for a three-dimensional system where the wavefunction for a system has the value 𝜓 at the point r with the coordinates (x, y, z), the probability of finding the particle in a miniscule volume 𝑑𝜏 = 𝑑𝑥𝑑𝑦𝑑𝑧 at the location r is proportional to |𝜓|2𝑑𝜏 [12].
To apply the Schrödinger equation in three dimensions, the Laplacian operator is added. The Laplacian operator is shown in equation (3) and the adapted Schrödinger equation in (4).
∇2= 𝜕 2 𝜕𝑥2+ 𝜕2 𝜕𝑦2+ 𝜕2 𝜕𝑧2 (3) − ħ 2 2𝑚∇ 2𝜓 + 𝑉𝜓 = 𝐸𝜓 (4)
A systematic way of extracting more information from the wavefunction, is to use the operator form of the Schrödinger equation: Ĥ𝜓 = 𝐸𝜓. In this equation, there is one new term in place of the previous expressions of kinetic and potential energy: Ĥ, the Hamiltonian operator. The purpose of the operator is to carry out a mathematical operation on the wavefunction. The Hamiltonian operator, or the electronic Hamiltonian, is defined in equation (5). The first term in the operator is the kinetic energy of the electrons Ne. The second term is the potential energy between each electron to each nucleus Nn, which is an attractive form of energy because of their opposite charges. From nucleus 𝐼 with the charge 𝑍𝐼𝑒, the electron 𝑖 is at a distance 𝑟𝐼𝑖. See Figure 1 for a schematic
6
Figure 1. A schematic description of the distance 𝑟𝐼𝑖 between the nuclei I=A and B and respective electron, as well as the distance 𝑟𝑖𝑗 between the electrons.
The third term takes the potential energy of repulsion between two electrons into consideration. The electrons are separated by 𝑟𝑖𝑗. The factor ½ is included in the third term to ensure that each
repulsion is only counted once [13].
Ĥ = − ħ 2 2𝑚𝑒 ∑ ∇𝑖2 𝑁𝑒 𝑖 − ∑ ∑ 𝑍𝐼 𝑒 2 4𝜋𝜀0𝑟𝐼𝑖 𝑁𝑛 𝐼 𝑁𝑒 𝑖 +1 2∑ 𝑒2 4𝜋𝜀0𝑟𝑖𝑗 𝑁𝑒 𝑖≠𝑗 (5)
In the first term of the Hamiltonian operator is the Laplace operator, ∇2. This operator, embedded into the Hamiltonian operator, makes sure that the system is studied in three dimensions. See equation (4) for the definition of the Laplace operator.
The term 𝑒2
4𝜋𝜀0 is a common reoccurrence in all computational chemistry and is henceforth indicated
as 𝑗0. With this change, the operator is as in equation (6):
Ĥ = − ħ 2 2𝑚𝑒 ∑ ∇𝑖2 𝑁𝑒 𝑖 − 𝑗0∑ ∑ 𝑍𝐼 𝑟𝐼𝑖 𝑁𝑛 𝐼 𝑁𝑒 𝑖 +1 2𝑗0∑ 1 𝑟𝑖𝑗 𝑁𝑒 𝑖≠𝑗 (6)
An interesting thing to note about the Schrödinger equation written in the operator form, it is that it is an eigenvalue equation. This means that the equation has the form (operator)(function) =
(constant factor)(same function). The eigenvalue equation can be formulated as 𝛺𝜓 = 𝜔𝜓, where 𝛺 represents the operator and 𝜔 is the constant factor. The 𝜔 is what is called the eigenvalue of the operator, which in the Schrödinger equation is the energy. The function 𝜓 is different for each eigenvalue 𝜔, and is the eigenfunction of the operator 𝛺. Considering the structure of the eigenvalue
7
equation on the Schrödinger equation, one can deduce that the eigenfunction is the wavefunction which corresponds to the energy E. So mathematically, to find the eigenvalues and the
eigenfunctions of the Hamiltonian operator is to solve the Schrödinger equation [12].
3.1.2 Molecular orbitals- Linear combinations of atomic orbitals
The Schrödinger equation can be solved analytically for hydrogenic atoms and the 𝐻2+molecule.
However, the wavefunctions are complicated types of functions, and analytical solutions cannot be obtained for other polyatomic systems. Therefore, an approximation can be made where a linear combination of atomic orbitals, LCAO, is constructed to apply to a molecular system. The general form of a LCAO for a polyatomic molecule is written in equation (7), where 𝜒𝑜represents an atomic
orbital and c is a coefficient [14].
𝜓 = ∑ 𝑐𝑜𝜒𝑜 𝑜
(7)
As previously mentioned in 3.1.1., 𝜓2 is proportional to the probability density of an electron at a location r.
3.1.3 Hartree-Fock theory
For a many electron molecule, the electronic wavefunction can be described as the function 𝛹(𝑟1, 𝑟2, … ). A common approximation is the orbital approximation, where the complete
wavefunction instead is written as 𝜓(𝑟1)𝜓(𝑟2) …. This approximation, however, is severe because
many details of how the wavefunction is dependent on the relative locations of the electrons are lost. In the following equations 𝜓(𝑟1)𝜓(𝑟2) … will be simplified as 𝜓(1)𝜓(2) …. For a given
molecule, suppose that electron 1 with the spin 𝛼 occupies the orbital 𝜓𝑎, and electron 2 has the
spin 𝛽 and also occupies the orbital 𝜓𝑎 and so on. This way, the many electron wavefunction can be
constructed as 𝛹 = 𝜓𝑎𝛼(1)𝜓𝑎 𝛽
(2) …. This function does not satisfy the Pauli principle, which states that the total wavefunction changes sign when the labels are exchanged between two identical fermions (fermions are particles with ½ spin: electrons, protons, and neutrons to name a few). Thus, to make sure that the function can satisfy the principle, it is modified to a sum of all the possible combinations: 𝛹 = 𝜓𝑎𝛼(1)𝜓𝑎 𝛽 (2) … 𝜓𝑧𝛽(𝑁𝑒) − 𝜓𝑎𝛼(2)𝜓𝑎 𝛽 (1) … 𝜓𝑧𝛽(𝑁𝑒) + ⋯ (8)
In the function presented above there are 𝑁𝑒! terms, and the total sum can be represented by the
Slater determinants: 𝛹 = 1 √𝑁𝑒 | | 𝜓𝑎𝛼(1) 𝜓𝑎𝛽(1) 𝜓𝑎𝛼 (2) 𝜓𝑎 𝛽 (2) … 𝜓𝑧𝛽(1) … 𝜓𝑧 𝛽 (2) ⋮ ⋮ 𝜓𝑎𝛼(𝑁𝑒) 𝜓𝑎 𝛽 (𝑁𝑒) . ⋮ … 𝜓𝑧𝛽(𝑁𝑒) | |
The Slater determinant is usually simplified by being written as the diagonal of the matrix:
𝛹 = (1/𝑁𝑒!)1/2|𝜓𝑎𝛼(1)𝜓𝑎 𝛽
(2) … 𝜓𝑧𝛽(𝑁𝑒) |
8
It is desired to find the optimum 𝜓𝑎 that will affect all other 𝜓 in the determinant, to create the best
form of 𝛹. The best form of 𝛹 is the one that achieves the lowest energy available. The aim of the equation (10) below, is to find 𝜓𝑎:
𝑓1𝜓𝑎(1) = 𝜀𝑎𝜓𝑎(1)
(10)
𝑓1 is the Fock operator which is an equation with the form 𝑓1= ℎ1+ 𝑉𝐶𝑜𝑢𝑙𝑜𝑚𝑏+ 𝑉𝐸𝑥𝑐ℎ𝑎𝑛𝑔𝑒. The first
term is the core Hamiltonian, meaning the one-electron Hamiltonian. The average Coulomb electron repulsion is represented by 𝑉𝐶𝑜𝑢𝑙𝑜𝑚𝑏, and the third term is the average correction concerning spin
correlation, 𝑉𝐸𝑥𝑐ℎ𝑎𝑛𝑔𝑒. The Coulomb repulsion contributes an integral, shown below in equation (11),
which represents the repulsion that electron (1) in the 𝜓𝑎 orbital experiences from electron (2) in the
orbital 𝜓𝑚. In the latter orbital, the probability density of the electron is 𝜓𝑚∗ 𝜓𝑚.
𝐽𝑚(1)𝜓𝑎(1) = 𝑗0∫ 𝜓𝑎(1)
1 𝑟12
𝜓𝑚∗ (2)𝜓𝑚(2)𝑑𝜏2
(11)
Since there are two electrons to each orbital, the total contribution of the Coulomb repulsion will be in the form:
𝑉𝐶𝑜𝑢𝑙𝑜𝑚𝑏𝜓𝑎(1) = 2 ∑ 𝐽𝑚(1)𝜓𝑎(1) 𝑚
(12)
In this form, one can see that the sum is including all the occupied orbitals, including 𝜓𝑎. However, it
is incorrect to count the orbital with m=a twice, since one electron in the orbital 𝜓𝑎 only interacts
with the other electron in this orbital and not with itself. This is corrected by adding the term
regarding average spin correlation, 𝑉𝐸𝑥𝑐ℎ𝑎𝑛𝑔𝑒. In this term, the avoidance between electrons of same
spin is taken into consideration, which then decreases the net Coulomb interaction. The form of the equation for the spin correlation is quite similar to the Coulomb repulsion and has the form:
𝐾𝑚(1)𝜓𝑎(1) = 𝑗0∫ 𝜓𝑚(1) 1
𝑟12𝜓𝑚∗ (2)𝜓𝑎(2)𝑑𝜏2
(13) The total contribution of this term, because there in one orbital only can reside two electrons of opposite spin, will be:
𝑉𝐸𝑥𝑐ℎ𝑎𝑛𝑔𝑒𝜓𝑎(1) = − ∑ 𝐾𝑚(1)𝜓𝑎(1) 𝑚
(14)
All these terms of the Fock operator can be gathered into one expression, which then describes the effect of the operator:
𝑓1𝜓𝑎(1) = ℎ1𝜓𝑎(1) + ∑{ 𝑚
2𝐽𝑚(1) − 𝐾𝑚(1)}𝜓𝑎(1)
9
As previously stated, to find 𝜓𝑎 the equation (10) needs to be solved. However, the equations (11)
and (13) proves it necessary to find all other occupied wavefunctions to be able to define the
operators J and K and through that find 𝜓𝑎. One approach to solve this is to initially guess the form of
all one-electron wavefunctions. The guessed forms are then used in the equations for Coulomb and exchange and are then used to solve the HF-equations. This process is then continued by using the wavefunctions given by one previous calculation in the following calculation, until the energies 𝜀𝑎
and wavefunctions 𝜓𝑎 are unchanged to certain criteria under several cycles. This type of process is
what is called a self-consistent field, SCF, and when used on the orbital approximation approach as described in this chapter, it is called Hartree-Fock self-consistent field, HF-SCF [13].
A problem when using the SCF procedure is the immense amount of numerical solutions of HF-equations. A way to work around this, and to make the technique more accessible, is to make modifications. This is done by applying the approach of molecular orbitals written as linear
combinations of atomic orbitals, LCAO, as described in chapter 3.1.2. Thus, as seen in equation (14), the coefficients of the LCAO that simulate the molecular orbitals are made by converting the HF-equations for molecular orbitals:
𝜓𝑚 = ∑ 𝑐𝑜𝑚𝜒𝑜 𝑁𝑏
𝑜=1
(14) In this equation, 𝑐𝑜𝑚 are the unknown coefficients, and 𝜒𝑜represents the atomic orbitals. When
using this linear combination, it will yield a set of equations for the coefficients called the Roothaan equations, which are commonly summarised as the matrix form:
𝐹𝑐 = 𝑆𝑐𝜀
(15) In this equation, F represents an 𝑁𝑏x𝑁𝑏 matrix as seen in equation (16). S is shown in equation (17)
and is an 𝑁𝑏x𝑁𝑏 matrix of integrals which are overlapping. c is also an 𝑁𝑏x𝑁𝑏 matrix, but one that
represents all the coefficients that are to be found, where each column in the matrix represents the coefficients for a specific 𝜓. The last term is a diagonal matrix, 𝜀, where the orbital energies are specified. 𝐹𝑜𝑜´ = ∫ 𝜒𝑜(1)𝑓1𝜒𝑜´(1)𝑑𝜏1 (16) 𝑆𝑜𝑜´ = ∫ 𝜒𝑜(1)𝜒𝑜´(1)𝑑𝜏1 (17) When writing the Roothaan equations as (𝐹 − 𝑆𝜀)𝑐 = 0, it can be deduced that they only have a solution when |𝐹 − 𝜀𝑆| = 0. In theory, the orbital energies of 𝜀 can be found by solving this equation for its roots, and then applying the energy values to the Roothaan equations to find the coefficients c. However, these equations are solved iteratively, where initially the values for c are guessed to solve the equation |𝐹 − 𝜀𝑆| = 0 for the orbital energies. The orbital energies are then used to solve
10
the Roothaan equations for c, which then are compared to the initial values of c. This cycle will then repeat until the values for c correlates [13].
A basis set, depicted as 𝑁𝑏 in this chapter, is used when constructing molecular orbitals by defining
which set of atomic orbitals that are to be used in the calculations [26]. The basis set is not
necessarily the same number as the number of atoms in a molecule, as different approaches can be made. The easiest approach is the minimal basis set, which uses one basis function per type of orbital included in the construction of the MO. With this approach, Ne for example would have five basis set functions: one for 1s and 2s respectively, and the three 2p-orbitals [13].
The HF-approach along with using basis set functions requires the calculation of a sizable number of integrals, making it very demanding regarding time as well as resources. To overcome this problem, several approximations can be made to the calculation of many-electron molecules. These
approximations are all different in nature and they each come with various sacrifices regarding accuracy [13].
3.2 Computational methods
3.2.1 First approach – ab initio methods
The so-called ab initio (“from the beginning”) methods originate from, or used as a first
approximation to, the Hartree-Fock theory. In Hartree-Fock calculations, however complicated they are, there is at least one aspect that can make the results less exact than one is expecting. Since the HF-method is based on orbital approximation and the average effect of the impact that one or more electrons have on the electron which is studied, the method does not take in consideration the fact that electrons are inclined to keep a certain distance to one another to minimize the repulsion between them. When the wavefunction is written as a product of functions 𝜓(𝑟1)𝜓(𝑟2) …., the
complexity of the movements of the electrons relative to one another is missing. There are two main approaches within the ab initio methods, and these are configuration interaction and many-body perturbation theory [13].
In configuration interaction, the aim is to improve the HF equation by including configuration state functions, CSF. A CSF is a Slater determinant which has been modified to include virtual orbitals 𝜓𝑣 as
well as the highest occupied molecular orbital, HOMO, 𝜓𝑢. By doing this, one can imagine an
electron being transferred from HOMO to a virtual orbital. The number of unoccupied virtual orbitals can be calculated as Nb- ½Ne, where Nb is the molecular orbitals and Ne is the number of electrons. Since electrons pair with opposite spin in the ground state of a molecule, only half of the number of electrons are deducted from the total number of orbitals. This ground state would then be written as:
𝛹
0= (1/𝑁
𝑒!)
1/2|𝜓
𝑎𝛼(1)𝜓
𝑎 𝛽(2)𝜓
𝑏𝛼(3)𝜓
𝑏𝛽(4) … 𝜓
𝑢𝛽(𝑁
𝑒)|
Where 𝜓𝑢 signifiesthe highest occupied molecular orbital, HOMO. From here one can visualise the
excited state where one electron has moved to a virtual orbital 𝜓𝑣:
𝛹
1= (1/𝑁
𝑒!)
1/2|𝜓
𝑎𝛼(1)𝜓
𝑎 𝛽(2)𝜓
𝑏𝛼(3)𝜓
𝑣𝛽(4) … 𝜓
𝑢𝛽(𝑁
𝑒)|
In principle, by constructing a linear combination of CSFs, they can be used to express the exact wavefunction of a system. For the calculations to be exact though, an infinite number of CSFs would be required [13].
In many-body perturbation, as the name reveals, perturbation theory is applied to the system of interacting nuclei and electrons within a molecule. In accordance with perturbation theory, the
11
Hamiltonian operator is expressed as a sum of a “model” Hamiltonian, Ĥ(0), and a perturbation Hamiltonian, Ĥ(1): Ĥ = Ĥ(0)+ Ĥ(1). The aim is to find the difference between the true energy and
the energy calculated by Hartree-Fock, and this is called the correlation energy. Because of its base in the HF-approach, the Fock-operators (equation 18) are commonly used as the “model” Hamiltonian, and the many-electron Hamiltonian, as seen in equation 6, as the perturbation Hamiltonian [13].
3.2.2 Second approach - density functional theory
Density functional theory is also based on the Schrödinger equation, but unlike previously named techniques, the aim of DFT is not to solve for the wavefunction. Instead, DFT can derive the distribution of electrons in a molecule immediately [11]. This is accomplished by solving for the electron probability density, ρ, through the energy functional E[ρ]. In mathematics, a functional is a function of another function, and in a molecular system this correlates to the fact that the energy of a molecule is a function of the density of its electrons, and that electron density is in itself a function of the positions of the electrons. From the structure of the Hamiltonian operator in equation 6, it can be deduced that the energy of a molecule can be described as kinetic energy, electron-nuclei
interaction, and electron-electron interaction. The contribution from the kinetic energy and the electron-nuclei interaction depend the density of the electron distribution, but because of electron exchange, the classic electron-electron interaction needs to be modified. The Hohenberg-Kohn equation shows how the contribution of electron exchange can be expressed in terms of electron density [13]:
𝐸[𝜌] = 𝐸𝐶𝑙𝑎𝑠𝑠𝑖𝑐𝑎𝑙[𝜌] + 𝐸𝑋𝐶[𝜌]
𝐸𝐶𝑙𝑎𝑠𝑠𝑖𝑐𝑎𝑙[𝜌] is the sum of the classical potential energy between electron-electron, the interaction
between nuclei-electron, and the contributions of kinetic energy. 𝐸𝑋𝐶[𝜌] is the energy of
exchange-correlation in which all the non-classical terms are accounted for, such as the effects of electron in terms of spin, as well as some corrections to the kinetic energy regarding electron-electron interactions.
The electron density ρ is calculated using an equation which uses the contribution from each electron in a molecule to express ρ:
𝜌(𝑟) = ∑ |𝜓𝑖(𝑟)|2 𝑁𝑒
𝑖=1
The term 𝜓𝑖 represents a so-called Kohn-Sham orbital, which is a solution to the Kohn-Sham
equation. This equation is based on the Schrödinger equation and is here presented as applied to a two-electron system with 𝜀𝑖 representing the Kohn-Sham orbital energies:
ℎ1𝜓𝑖(1) + 𝑗0∫
𝜌(2) 𝑟12
𝑑𝜏2𝜓𝑖(1) + 𝑉𝑋𝐶(1)𝜓𝑖(1) = 𝜀𝑖𝜓𝑖(1)
The first term in this equation is the core term, and the second term represents the classical interaction between the two electrons. The third term is what is called an exchange-correlation potential, which accounts for the effects of exchange. To solve this exchange-correlation term of the Kohn-Sham equation, it can be expressed as a functional derivative of the exchange-correlation energy:
12 𝑉𝑋𝐶(𝑟) =
𝛿𝐸𝑋𝐶[𝜌]
𝛿𝜌
These equations are solved iteratively and self-consistently. The first step is usually to use a superposition of atomic electron probability densities to guess the electron density. The following step is to calculate the exchange-correlation potential by approximating a form of dependence of the exchange-correlation energy using an assumption and then assessing the functional derivative. The next course of action is to find a set of Kohn-Sham orbitals by solving the Kohn-Sham equations, this set of orbitals is then used to obtain a better approximation of the electron probability density, and the calculation process is then repeated until the density remains constant [13].
DFT is a method that is more computationally efficient in comparison to HF and is generally the preferred method in which to calculate molecular structures [13].
4. Method
4.1 Computer and programs
Calculations were carried out using the supercomputer Tetralith on the national supercomputer, NSC, which is located in Linköping and is a part of Linköping University and the Swedish National Infrastructure for Computing, SNIC. Tetralith contains a total of 1908 compute nodes. Each node contains two CUPs with 16 CPU cores each [15]. To access the supercomputer, the remote desktop ThinLinc was used. The program Gaussian was used together with the program GaussView 6.1.1. to build molecules and to set up calculation parameters. The command sgausbatch was used when submitting jobs for calculation on Tetralith, whereby allocating 32 processors and 16 GB memory. The chosen computational level of theory for the calculations, both optimisations and single-point energies, was DFT. The settings were B3LYP (Becke three-parameter Lee-Yang-Par) exchange-correlation functional, which is a gradient-corrected HF/DFT hybrid functional [16]. The basis set 6-31G(d,p) was used, which is a split valence basis set with the polarisation function (d,p). Split valence means that the valence shell of each atom is split in two, where one inner part is made up of three functions and an outer part consisting of one function. The core orbital for each atom is
characterized by one basis function which is made up by six functions. The polarisation (d,p) is used on atoms heavier than He [17]. The empirical dispersion correction used for all calculations was set to GD3.
To recover files such as output-files or pictures from Tetralith and ThinLinc, the program WinSCP was used.
4.2 ZnO-particles: size and structure
Initially, the appropriate size of the ZnO particle was investigated with the aim to construct a particle that would yield a low energy of formation in comparison to the energy of the separate atoms, which can be called the cohesive energy. This energy was calculated using the formula written below in formula (18):
∆𝐸𝑐𝑜ℎ= 𝐸((𝑍𝑛𝑂)𝑛) − (𝑛𝐸(𝑍𝑛) + 𝑛𝐸(𝑂)) (18)
A negative value indicates bonding. The particles were constructed using GaussView. The Zn and O atoms were uncharged atoms, and by connecting n number of Zn-atoms with n number of O-atoms, structures of different sizes were formed while ensuring as few dangling bonds as possible.
13
The particles constructed were of the size (ZnO)4 to (ZnO)20, and several different types of structures were made. A selected number of particles were created in two versions, versions 1 and version 2, to be able to compare the energy of formation of the structures. When the structures had been
created, a geometry optimisation was conducted. The geometry optimisation was performed to find the geometrical structure that had the lowest energy regarding the attracting and repelling forces within the ZnO-cluster.
The energies of the respective ZnO-clusters were converted into the unit kJ/mol instead of Hartree (1 Hartree= 2625,5 kJ/mol), and the difference in energy ∆𝐸 of each system was calculated using the formula (18). ∆𝐸 for each particle was then divided with the number n of (ZnO)n and the energy was thereafter plotted in a diagram on the y-axis with n on the x-axis.
4.3 ZnO-particles: oxygen vacancy
The particle with the lowest ∆𝐸𝑐𝑜ℎ was then used in further analysis. One oxygen atom in the
structure was removed, at different locations in the particle, creating new clusters ZnnOn-1 with different structures. Geometrical optimisations and energy calculations were performed on the new structures as well to identify from which cluster an O atom could be removed while requiring the least amount of energy. The same optimisation and energy set up was used as in prior calculations. The formula to calculate the difference in energy between the complete (ZnO)n structure with the particle with one O atom removed, i.e. the vacancy formation energy, was:
∆𝐸𝑣𝑎𝑐 = (𝐸(𝑍𝑛𝑛𝑂𝑛−1) + 𝐸(𝑂)) − 𝐸(𝑍𝑛𝑂)𝑛 (19)
4.4 CO
2bonding to ZnO
After identifying which clusters would be more beneficial regarding energy, the ability to adsorb CO2 molecules at the surface of said clusters was investigated. This was done by introducing a CO2 molecule near the particles and performing a geometry optimisation to find the Zn atom in the structure to which one oxygen from CO2 would coordinate by physisorption most favourably. This was investigated on both (ZnO)n and ZnnOn-1, again with the purpose to compare the species. Before the optimisation was performed on the (ZnO)n particle, it was necessary to find a location near the cluster where to strategically place the CO2 molecule in an initial position. In GaussView, the molecular orbitals for the cluster was studied visually, and the CO2 was placed close to the location where the highest occupied molecular orbital, HOMO, was found. For the ZnnOn-1 particle, the CO2 was placed close to the oxygen vacancy.
After the geometrical optimisation, the energy was calculated. The settings used for the calculations were as described in chapter 4.1.
To be able to compare the clusters, the difference in energy of the system upon adsorption of CO2 at a cluster with an O vacancy or without an O vacancy was calculated as given in the formulas (20) and (21).
∆𝐸𝑎𝑑𝑠𝑣𝑎𝑐 = 𝐸(𝑍𝑛𝑛𝑂𝑛−1+ 𝐶𝑂2) − (𝐸(𝐶𝑂2) + 𝐸 (𝑍𝑛𝑛𝑂𝑛−1)) (20)
14
4.5 Activation energy
Next the ability of the ZnO particle to chemically react with a CO2 molecule was investigated. To calculate the activation energy required for CO2 to react with a particle of ZnO, the particle that showed the most preferable regarding energy in previous experiments was used, meaning the cluster that yielded the lowest ∆𝐸𝑐𝑜ℎ and the lowest ∆𝐸𝑣𝑎𝑐. CO2 was placed in a relatively close vicinity to
the cluster, while freezing the coordinates for one Zn-atom at a strategic place in the cluster and one O-atom in CO2. An optimisation was performed using to the same settings as previously mentioned in 4.1. After the initial optimisation of CO2 at a distance from the cluster, the structures were optimised yet again, using a new position of CO2 relative to the cluster as the new starting point. This approach was applied until one oxygen from CO2 was detached and instead formed chemical bond to the cluster. When this occurred, energy calculations were performed on the reactant, the transition state, and the product. These energies were then used to calculate ΔE as described in formulas (22) for the reaction to the right, and (23) for the reaction to the left. A diagram was created where ΔE was visualised as a function of the bond length of C-O (the O-atom which is approaching Zn), which was chosen as the reaction coordinate (x-axis), from which the activation energy could be identified.
∆𝐸𝑎ⴕ= 𝐸(𝑡𝑟𝑎𝑛𝑠𝑖𝑡𝑖𝑜𝑛 𝑠𝑡𝑎𝑡𝑒) − 𝐸(𝑟𝑒𝑎𝑐𝑡𝑎𝑛𝑡) (22)
∆𝐸𝑎ⴕ= 𝐸(𝑡𝑟𝑎𝑛𝑠𝑖𝑡𝑖𝑜𝑛 𝑠𝑡𝑎𝑡𝑒) − 𝐸(𝑝𝑟𝑜𝑑𝑢𝑐𝑡) (23)
5. Results
5.1 Size and structure of ZnO
The structures created were of the size n=4-16 and n=18-21, whereas n=17 was excluded because of its repeatedly failing geometry optimisation. A few clusters, n= 4, 9, 13, 15, 19 and 20, were created in two different versions to see whether there would be a difference in energy depending on the structure. A selected few of the optimised structures, n=10, 14, 19 and 20 version 2, are shown in Figure 2, where the red atoms represent O, and the blue represent Zn. All clusters can be seen in the appendix, Figure 1-12.
Figure 2: Optimised structures of (ZnO)n. Left: Top (ZnO)10, bottom (ZnO)19 version 2. Right: Top (ZnO)14, bottom (ZnO)20
15
Energy calculations were performed for the geometrically optimised clusters, the corresponding energy data for each cluster of (ZnO)n in the unit Hartree as well as converted to kJ/mol is presented in Table 1. The energies for Zn and O, along with the calculation of ΔEcoh, using formula (18), for each cluster and ΔEcoh/n is also presented in this table, in the third and fourth column labelled ΔEcoh and ΔEcoh/n.
Table 1: Energy (Hartree and kJ/mol) for each cluster (ZnO)n, and the calculated values for ΔEcoh and ΔEcoh/n.
ZnO-particle size (n)
E (Hartree) E (kJ/mol) ΔEcoh (kJ/mol) ΔEcoh/n (kJ/mol)
Zn -1779,09 -4 671 001,24 O -75,06 -197 071,66 4 version 1 -7417,53 -19 474 715,83 -2424,23 -606,06 4 version 2 -7417,47 -19 474 574,25 -2282,66 -570,66 5 -9271,874 -24 343 290,44 -2925,95 -585,19 6 -11 126,38 -29 212 316,00 -3878,60 -646,43 7 -12 980,80 -34 081 097,41 -4587,11 -655,30 8 -14 835,21 -38 949 840,49 -5257,29 -657,16 9 version 1 -16 689,62 -43 818 608,11 -5952,01 -661,33 9 version 2 -16 689,71 -43 818 834,29 -6178,19 -686,47 10 -18 544,13 -48 687 619,51 -6890,52 -689,05 11 -20 398,55 -53 556 383,64 -7581,75 -689,25 12 -22 253,01 -58 425 280,65 -8405,86 -700,49 13 version 1 -24 107,43 -63 294 059,40 -9111,70 -700,90 13 version 2 -24 107,46 -63 294 143,41 -9195,72 -707,36 14 -25 961,89 -68 162 937,23 -9916,64 -708,33 15 version 1 -27 816,30 -73 031 696,62 -10 603,13 -706,88 15 version 2 -27 816,37 -73 031 884,01 -10 790,52 -719,37 16 -29 670,76 -77 900 590,46 -11 424,07 -714,00 18 -33 379,64 -87 638 245,34 -12 933,14 -718,51 19 version 1 -35 234,05 -92 507 005,32 -13 620,23 -716,85 19 version 2 -35 234,15 -92 507 255,17 -13 870,08 -730,00 20 version 1 -37 088,52 -97 375 900,32 -14 442,33 -722,12 20 version 2 -37 088,62 -97 376 180,30 -14 722,31 -736,12 21 -38 943,03 -102 244 922,90 -15 392,04 -732,95
ΔEcoh/n and n from Table 1 was plotted in a diagram, see Figure 3, to visualise which value of n that provides the lowest ΔEcoh. A trendline is also presented in the diagram, as well as the third degree polynomial which describes the line and thus the correlation between x=n and y=ΔE(kJ/mol)/n, which is 𝑦 = −0,0537𝑥3+ 2,61𝑥2− 45,014𝑥 − 442,54. The diagram shows a fairly monotonous relation between ΔEcoh/n and n, where ΔEcoh decreases with larger n.
16
Figure 3: A diagram of (ZnO)n-particles of different size n, and their respective calculated value of ΔE/n. A trendline is customised to the data in the diagram and a third degree polynomial is shown which describes the line.
According to Table 1 and Figure 3, the structure that provides the lowest ΔEcoh is the second version of (ZnO)20, the structure of which is seen at bottom right in Figure 2. Both version 1 and 2 of (ZnO)20 was then used to research O vacancies.
5.2 Oxygen vacancy
An O-atom was eliminated at different locations in the structure for version 1 of (ZnO)20, creating three different structures of Zn20O19. In version 1:1 an O-atom was removed at the end of the structure, while in 1:2 and 1:3 the atoms removed were in the middle of the structure. For version 2 of (ZnO)20, the oxygen was removed from the centre of the largest surface. The Zn20O19-structures, four in total, were geometrically optimised, and an energy calculation was performed. The data for the energies and ΔEvac calculated using formula (19) are shown in Table 2.
Table 2: Energy (Hartree and kJ/mol) and the calculated values for ΔEvac and ΔEvac/n. The three structures of Zn20O19 for
version 1 of (ZnO)20, as well as Zn20O19 for version 2.
Zn20O19
-structure
E (Hartree) E (kJ/mol) ΔEvac (kJ/mol)
Version 1:1 -37 013,22 -97 178 216,09 612,57
Version 1:2 -37 013,21 -97 178 171,84 656,82
Version 1:3 -37 013,20 -97 178 165,63 663,02
Version 2:1 -37 013,30 -97 178 418,03 410,62
Table 2 shows that version 2:1 of Zn20O19 has the lowest calculated ΔEvac of all. Of the three versions of structure 1, the version 1:1 has the lowest ΔEvac. The two versions, 1:1 and 2:1, are pictured in Figure 4 where the difference in structure and the different placements of oxygen vacancy can be seen. y = -0,0537x3+ 2,61x2- 45,014x - 442,54 -750 -700 -650 -600 -550 -500 0 5 10 15 20 25 Δ E(kJ /m o l)/ n n
17
Figure 4: Oxygen-vacancies in (ZnO)20-particles, forming Zn20O19. To the left: version 2:1 from above with vacancy in the
middle of the top layer of the structure as indicated by the arrow. To the right: Version 1:1 from the side, with oxygen vacancy to the far left as indicated by the arrow.
Because of their lower ΔEvac, the versions 1:1 and 2:1 were used in the following calculations where ΔEads of the bonding of CO2 to the clusters was investigated, along with the two original versions of (ZnO)20 without O vacancy.
5.3 CO
2adsorption
The highest occupied molecular orbital, HOMO for version 2 of the (ZnO)20 cluster is presented in figure 5. The arrow in the figure indicates where CO2 was placed ahead of optimising the structure. Similarly, the HOMO was investigated for version 2 of (ZnO)20 and CO2 was thereafter placed strategically. For the structures with an O vacancy, versions 1:1 and 2:1 of Zn20O19, the CO2 was placed near the vacancy, see Figure 4 for clarification as to where these are located. In total, four structures with CO2 were optimised.
Figure 5: A picture of the highest occupied molecular orbital, HOMO, of (ZnO)20. The arrow indicates the vicinity where CO2
was added prior to optimisation.
The results of the energy calculation after optimisation of the four structures are presented in Table 3. This table also contains the calculated value of ΔEads, which is presented in the fourth column. When comparing ΔEads depending on whether the structure has an O-vacancy or not, the first version of (ZnO)20 gives quite similar values of ΔEads for both, while the second version differs significantly in ΔEads.
18
Table 3: Energy in Hartree and kJ/mol and calculated ΔEads for versions 1 and 2 with and without O vacancy.
Structure E (Hartree) E (kJ/mol) ΔEads
CO2 -188,58 -495 119,75 - Version1: (ZnO)20+CO2 -37 277,12 -97 871 081,05 -60,97 Version 1:1: Zn20O19+CO2: -37 201,82 -97 673 388,82 -52,98 Version 2: (ZnO)20+CO2 -37 277,26 -97 871 441,32 -141,27 Version 2:1: Zn20O19 + CO2 -37 201,88 -97 673 528,92 8,87
Table 3 shows that version 2, which is the second version of (ZnO)20 without an O vacancy, is the most energy sufficient when adsorbing CO2. However, further research will be performed on the version 2:1, with an O vacancy, because the presence of an O vacancy is probably necessary for the catalytic action of ZnO where in a reaction step in the methanol synthesis an oxygen is removed from CO2 [18].
5.4 Conversion of CO
2to CO
The activation energy for the version 2:1 Zn20O19 cluster was investigated by identifying the reactant, the transition state, and the product of the reaction Zn20O19+CO2 →(ZnO)20+CO. The coordinates of a Zn atom in direct vicinity to the oxygen vacancy along with one O atom of CO2 were frozen during the geometry optimisation. See Figure 6, where the arrows indicate the Zn and O atoms used for freezing coordinates. The coordinates were used in pairs which were tested one at a time, these pairs were the atoms labelled as 41-25, 42-25, 42-15 and 41-15. The coordinates 42-15 proved to be the most successful in optimisation and used thereafter. Geometry optimisations were performed on a series of clusters with said coordinates frozen, which led to CO2 slowly approaching the cluster Zn20O19.
Figure 6: A picture of Zn20O19 and CO2 where all atoms are labelled. Zn atoms are blue and O atoms red. Zn with the
19
The initial optimised structure was used as the reactant structure. The optimised structure where one O atom, number 42, from CO2 had created a bond with one Zn atom, number 15, in close vicinity to the O vacancy in the cluster was used as the transition state. The product structure was the structure when an O atom had released, and CO was formed. Figure 7 shows the optimised structures for the transition state (top two pictures) and the product (bottom two pictures).
Figure 7: Left picture shows the structures from the side, and the right picture is the structures from above. Top: Zn20O19
with CO2 attached, this was used as the transition state. Bottom: (ZnO)20 with CO attached, this was used as the product.
The calculated energies are shown in Table 4. The calculated values for ΔE are also presented here, where the ΔE for the transition state is calculated for the reaction to the right, i.e. the formula (22). ΔE for the product was calculated using formula (23) for the reaction to the left. The ΔE for the reactant was set to 0, since it is the difference in energy from this structure that is desired.
Table 4: The calculated energies for the reactant, transition state and product in the reaction Zn20O19+CO2 →(ZnO)20+CO
Structure Length C-O (Å) E (Hartree) E (kJ/mol) E-E(reactant)
Reactant 1,167 -37 201,90 -97 673 590,39 0
Transition state 1,186 -37 201,85 -97 673 463,10 127,28
20
When plotting these values in a diagram, see Figure 7, the activation energy Ea can be read as the top of the curve, which gives the Ea= 127 kJ/mol.
Figure 7: A diagram of the activation energy required for the reaction Zn20O19+CO2 →(ZnO)20+CO with ΔE (kJ/mol) on the
y-axis and the distance between C-O (Å) as the reaction coordinate.
6. Conclusion and discussion
The aim of this study was to investigate by quantum-chemical computations the catalyst component ZnO used in the synthesis of methanol from CO2 and H2 using heterogeneous catalysis. In this study, the most energy-sufficient sizes and structures of ZnO-particles were therefore studied as well as the formation of an O vacancy on the surface of the cluster. The adsorption of CO2 was studied on ZnO with and without O vacancies, and lastly the activation energy required when forming CO from CO2 was calculated.
The optimisation of different sizes of (ZnO)n proved, within n=4-21, that n=20 was the better size in that it gave the lowest energy of formation. While comparing structures, one can see that the energies of versions 1 and versions 2, where the version 2 was structurally more similar to a lattice structure, differs in all cases where this was applied. This proves that the structure, while perhaps not as significant as the number n, is still relevant when optimising these particles. One can discuss whether the true most optimal n lies above n=21, but in the aspect of the limited time of this project, this matter is not investigated further since it would demand a lot more time.
Regarding oxygen vacancy, the structure of the ZnO particle proved even more important. There was a significant difference in energy of the optimised structures depending on the structure. The
energies of version 1 vary slightly depending on where the O vacancy was located, but they are still quite similar. The energy of version 2 however, which better represents the surface of a ZnO particle, differ a lot from version 1. The optimised structure of version 2 with an O vacancy (2:1) proved to give a total amount of energy which was more than 200 kJ/mol less than the three structures of version 1 (1:1, 1:2 and 1:3). -350 -300 -250 -200 -150 -100 -50 0 50 100 150 1,165 1,17 1,175 1,18 1,185 1,19 1,195 Δ E (k J/m o l) Distance C-O (Å)
Activation energy
21
When comparing the difference in energy when adsorbing CO2 to (ZnO)20 and Zn20O19, the results vary depending on which version is studied. The first version shows a relatively small difference between the species, and it seems that whether it has an O vacancy or not might be less important with this structure. In version two however, the difference in energy of adsorption is large, where the structure without an O vacancy gives a much lower energy, meaning that the cluster (ZnO)20 with adsorbed CO2, has a lower energy than the two species separately. The reason for the difference between versions 1 and 2 might be because of the more stable structure of the second version, which is more similar to a lattice structure than version 1. I had expected the outcome to be different for version 2:1 and 2, thinking that the vacant structure (2:1) would more readily bond with a new O atom from the desire to become more structurally stable again, which then would make it a better structure to use as a catalyst. The result that was however the opposite, which makes it even more interesting and necessary to compare different structures. It becomes clear that structure plays such an important part in chemical reactions.
The calculated activation energy for Zn20O19 + CO2 was 127 kJ/mol. However, with a dissociation energy of 1072 kJ/mol for CO2, this seems like quite a small amount of energy. In the papers that have been used for research in this study, I have not managed to find a calculated value of the activation energy for this reaction, so it has unfortunately not been possible to find a value to compare with. It would also have been preferable if activation energies were calculated for all four structures, namely version 1 with and without an O vacancy and version 2 with and without an O vacancy, to be able to compare and come to a more well-founded conclusion. However, there simply was not enough time to complete this.
Regarding the following reaction steps in the methanol synthesis, it is worth noting that ZnO with an O vacancy as a catalyst is a part of the reaction, but it cannot be consumed or changed in any way. This requires a recreation of the O vacancy before the reaction is stopped. The oxygen which bonds into the vacancy and subsequently is removed from CO2 forming CO, will need to create a new bond with another species to recreate the O vacancy. Possible reaction paths are that two O atoms in adjacent vacancies form O2(g), or that one O forms H2O with H2(g) which is a reactant in the CH3OH formation reaction. Still, the formation of CO is just a step in the reaction of forming methanol from CO2, and the molecule will need to be reduced further with H2 using the interaction with the catalysts Cu and ZnO.
The aim of this study was to research the most energy sufficient way to synthesise methanol from CO2 in regard to the catalyst component ZnO in nanoparticle form. This study concludes that the better size is (ZnO)20 with a large, flat surface, as it has a slight resemblance to a lattice structure. The same (ZnO)20 structure is the one that yields the least amount of energy when removing an O atom from the cluster. When adsorbing CO2 however, the structure without the O vacancy proved to be more favourable to engage in physisorption with CO2. The activation energy was investigated for the structure that proved better with an O vacancy, and was calculated to be 127 kJ/mol. To properly gauge the activation energy, and to draw more conclusions as to the better structure, size, and O vacancy of ZnO, further investigation is required.
Using computational methods to develop and optimise the synthesis of methanol through
heterogeneous catalysis of CO2 is, while quite time consuming, a very clever approach to the problem at hand regarding various aspects. One aspect is economical, where a full scale practical optimisation of this type of method would most likely be expensive, the computational approach offers a
theoretical optimisation of the method which can help to narrow down the possible solutions and to identify which parameters of a particular problem might be more suitable to test practically. The
22
same applies to the aspect of time, while the computational methods take a fair amount of time, to optimise the method in a practical setting would be even more time consuming.
A slightly negative aspect which has been bothersome in this project has been the failing
calculations. It has been difficult to predict whether the calculations, optimisations predominantly, would be successful and being carried out to completion or not. This led to many re-calculations, and many times it caused me to build new structures to calculate. However, this probably stems from my personal inexperience with, and lack of knowledge regarding, computational chemistry and the programming used in this project. Regardless, computational chemistry is an invaluable tool in this field of research, and without which I would not have been able to reach the conclusions of this project. Still, to the future of methanol synthesis at large, further time and effort is required to fully optimise the method and its parameters.
7. Acknowledgements
A massive thank you to examiner Lars Ojamäe, and tutor Karl Rönnby for their patience, guidance, and expertise, without whom this thesis never would have been written. I would also like to acknowledge the resources made available by SNIC and NSC.
23
References
[1] Naturvårdsverket. Utsläpp i siffror: koldioxid [Internet]. Stockholm: Naturvårdsverket;2017. Available from: https://utslappisiffror.naturvardsverket.se/Amnen/Vaxthusgaser/Koldioxid/
[2] W. Wang, S. Wang, X. Ma, J. Gong. Recent advances in catalytic hydrogenation of carbon dioxide. Chem. Soc. Rev., 2011, 40, 3703-3727.
[3] NASA. Global climate change; carbon dioxide [Internet]. USA: California Institute of Technology; 2020. Available from: https://climate.nasa.gov/vital-signs/carbon-dioxide/. (Measurement is from march, 2020.)
[4] E. Byers, M. Gidden, D. Leclère, J. Balkovic, P. Burek, K. Ebi, P. Greve, D. Grey, P. Havlik, A. Hillers, et al. Global exposure and vulnerability to multi-sector development and climate change hotspots. 2018, Environ. Res. Lett. 13 055012. Available from: https://iopscience.iop.org/article/10.1088/1748-9326/aabf45
[5] Y. Sun, C. Huang, L. Chen, Y. Zhang, M. Fu, J. Wu, D. Ye. Active site study of Cu / Plate ZnO model catalysts for CO2 hydrogenation to methanol under the real reaction conditions. Journal of CO2
utilization 37; 2020 (55-64).
[6] P. Atkins, J. de Paula. Atkins’ Physical Chemistry, 10th edition. United Kingdom: Oxford University
press; 2014. Table 2C.5 p 971-975
[7] S. Dang, H. Yang, P. Gao, H. Wang, X. Li, W. Wei, Y. Sun. A review of research progress on
heterogeneous catalysts for methanol synthesis from CO2 hydrogenation. Catalysis Today 330: 2019
(61-75).
[8] P. Atkins, J. de Paula. Atkins’ Physical Chemistry, 10th edition. United Kingdom: Oxford University
press; 2014. Chapter 20: Chemical kinetics.
[9] P. Atkins, J. de Paula. Atkins’ Physical Chemistry, 10th edition. United Kingdom: Oxford University
press; 2014. Chapter 22: Processes on solid surfaces.
[10] G. A. Somorjai. Introduction to surface chemistry and catalysis, 1st edition. USA: John Wiley &
sons, Inc; 1994. p 6-8.
[11] E. Lewars. Computational chemistry – the introduction to the theory and applications of
molecular and quantum mechanics. USA: Kluwer Academic Publishers; 2003. Chapter 1: An outline of what computational chemistry is all about.
[12] P. Atkins, J. de Paula. Atkins’ Physical Chemistry, 10th edition. United Kingdom: Oxford University
press; 2014. Chapter 7: Introduction to quantum theory.
[13] P. Atkins, J. de Paula, R. Friedman. Quanta, Matter and Change – a molecular approach to physical chemistry. United Kingdom: Oxford University press; 2009. Chapter 6: Computational chemistry.
[14] P. Atkins, J. de Paula. Atkins’ Physical Chemistry, 10th edition. United Kingdom: Oxford University
press; 2014. Chapter 10: Molecular structure.
[15] National Supercomputer Centre, NSC. Systems: Tetralith [Internet]. Linköping; 2020. Available from: https://www.nsc.liu.se/systems/tetralith/
24
[16] E. Lewars. Computational chemistry – the introduction to the theory and applications of molecular and quantum mechanics. USA: Kluwer Academic Publishers; 2003. Chapter 7: Density functional calculations.
[17] E. Lewars. Computational chemistry – the introduction to the theory and applications of molecular and quantum mechanics. USA: Kluwer Academic Publishers; 2003. Chapter 5: Ab initio calculations.
[18] J. Kobmann, G. Romüller, C. Hättig. Prediction of vibrational frequencies od possible
intermediates and side products of the methanol synthesis on ZnO (000) by ab initio calculations. J. Chem. Phys. 136, 034706 (2012).
25
Appendix
Figure 1: (ZnO)4 version 1 to the left and version 2 to the right.
Figure 2: (ZnO)5 to the left and (ZnO)6 to the right.
Figure 3: (ZnO)7 the left and (ZnO)8 to the right.
26
Figure 5: (ZnO)10 to the left and (ZnO)11 to the right
Figure 6: (ZnO)12 to the left and (ZnO)14 to the right
Figure 7: (ZnO)13 version 1 to the left and version 2 to the right.
27
Figure 9: (ZnO)16 to the left and (ZnO)18 to the right
Figure 10: (ZnO)19 version 1 to the left and version 2 to the right.
Figure 11: (ZnO)20 version 1 to the left and version 2 to the right.
