Multiscale Modeling of Electrolytes for Energy Storage and Conversion A Dissertation Presented for the Doctor of Philosophy Degree The University of Tennessee, Knoxville Fatemeh Sepehr December 2016 ii DEDICATION This dissertation is dedicated to my husband, my parents, and my sibling who have supported and encouraged me throughout my life and academic career. iii ACKNOWLEDGEMENTS I would like to express my sincere gratitude to my advisor, Dr. Stephen J. Paddison, for his tremendous support throughout my Ph.D. studies at the University of Tennessee. His guidance played an important role in my research as well as my professional growth. I sincerely thank the efforts of knowledgeable faculty whom I learned from during my academic endeavor. I would also like to extend my thanks to Dr. Chen Wang, Dr. Jeffrey Clark, Dr. Hongjuan Zhu, Dr. Jun He, Dr. Hongjun Liu, Dr. Dongsheng Wu, Mr. Vincent Hall, and Mr. Xubo Luo for their help. The financial support from the National Science Foundation (NSF) and the TNSCORE program is greatly appreciated and valued. iv ABSTRACT Fuel cells, redox flow batteries, and secondary ion batteries are under active investigation to fulfill the requirements of efficient and sustainable energy storage and conversion technologies. The discovery of high-performance stable electrolytes that are relatively cheap and versatile is crucial to the commercialization of these electrochemical devices and necessitates a comprehensive understanding of the materials (i.e., from the atomistic to continuum levels). This dissertation is on multiscale modeling and simulations of several electrolytes under consideration in vanadium redox flow batteries (VRFBs), alkaline fuel cells (AFCs), or secondary magnesium batteries. The hydrated structure and associated solvation Gibbs energies were determined for vanadium cations at four oxidation states using first principles based electronic structure calculations and quasi-chemical theory. The effect of sulfuric and trifluoromethanesulfonic (triflic) acids on the local structures of the hydrated vanadium cations was explored by employing hybrid density functional theory (DFT) in conjunction with a continuum solvation model. The results indicate that the hydration structure of vanadium at the oxidation state of three was the most perturbed in acidic media and that the oxo-group of vanadate may be protonated by either acid. Dissipative particle dynamics (DPD) simulations were undertaken and a novel parametrization method developed to study the meso-scale behavior and hydrated morphology of proton and anion exchange membranes (PEMs and AEMs). The results for hydrated Nafion with the absorbed vanadium species indicate that the hydrated cations increase the aggregation of the sulfonate groups and that this effect is a function of both the cation charge and the degree of hydration. It was also observed that the simulated morphology and the size of the water domains in AEMs based on triblock copolymers depend on the hydration level and the anion type. v Vibrational frequencies of model clusters were calculated using DFT for two ionic liquid-based electrolytes and their experimental IR and Raman spectra were assigned. Moreover, the hydration of ten acids commonly used in electrolytes was investigated using DFT in order to provide a fundamental understanding of their hydrated acidity and proton transfer. vi TABLE OF CONTENTS Chapter 1 Introduction .............................................................................................................1 1.1 Redox Flow Batteries ......................................................................................................1 1.1.1 Vanadium Redox Flow Batteries (VRFBs) .............................................................2 1.1.2 Electrolytes for VRFBs ...........................................................................................4 1.2 Fuel Cells .........................................................................................................................6 1.2.1 Proton Exchange Membrane Fuel Cells (PEMFCs) ..............................................6 1.2.2 Electrolytes for PEMFCs .......................................................................................8 1.2.3 Alkaline Fuel Cells (AFCs) ..................................................................................12 1.2.4 Electrolytes for Anion Exchange Membrane Fuel Cells ......................................13 1.3 Secondary Magnesium Batteries ...................................................................................16 1.3.1 Electrolytes Based on Ionic Liquids .....................................................................18 1.4 Research Overview ........................................................................................................19 Chapter 2 Computational Methods .......................................................................................22 2.1 Multiscale Modeling ......................................................................................................22 2.2 Quantum Chemical Calculations ...................................................................................24 2.2.1 Time-Independent Schrödinger Equation .............................................................25 2.2.2 Hartree-Fock Method ...........................................................................................28 2.2.3 Density Functional Theory ...................................................................................31 2.2.4 Basis Sets ..............................................................................................................35 2.2.5 Continuum Solvation Models ...............................................................................37 2.2.6 Vibrational Frequencies .......................................................................................40 vii 2.3 Quasi-Chemical Approach ............................................................................................42 2.3.1 Potential Distribution Theory ...............................................................................43 2.3.2 Basic Quasi-Chemical Formula ...........................................................................49 2.3.3 Chemical Equilibrium and Quasi-Chemical Description ....................................51 2.4 Dissipative Particle Dynamics .......................................................................................53 2.4.1 Determining the Interaction Parameters..............................................................55 Chapter 3 PFSA Electrolytes with Vn+ ..................................................................................59 3.1 Hydration and Thermodynamics of the Vn+ Cations .....................................................60 3.1.1 Hydration Structures ............................................................................................60 3.1.2 Gibbs Energies of Hydration ................................................................................63 3.2 Effect of Sulfuric and Triflic Acids on the Hydration of Vn+ ........................................67 3.2.1 The V2+ Cation ......................................................................................................68 3.2.2 The V3+ Cation ......................................................................................................70 3.2.3 The VO2+ Cation ...................................................................................................72 3.2.4 The VO2 + Cation ...................................................................................................74 3.2.5 Partial Atomic Charges ........................................................................................76 3.3 Approach to Parameterize Dissipative Particle Dynamics ............................................79 3.4 Morphology of the Membrane Electrolyte ....................................................................83 3.4.1 Coarse-Graining of Particles ...............................................................................84 3.4.2 Conservative Interaction Parameters ...................................................................85 3.4.3 Bonded Interaction Parameters ...........................................................................87 3.4.4 Simulations ...........................................................................................................89 Chapter 4 Anion Exchange Electrolytes ................................................................................94 viii 4.1 Coarse-Grained Model ..................................................................................................94 4.1.1 Coarse-Graining of the Particles .........................................................................94 4.1.2 Conservative Interaction Parameters ...................................................................96 4.1.3 Bonded Interaction Parameters ...........................................................................98 4.1.4 Simulations ...........................................................................................................99 4.2 Morphology of SEBS Triblock Copolymer ................................................................100 4.3 Morphology of SEBS-Based Electrolytes ...................................................................102 4.3.1 The Effect of Hydration ......................................................................................104 4.3.2 The Effect of Anion (OH− versus Cl−) ................................................................109 Chapter 5 SDAPP Electrolytes .............................................................................................116 5.1 Coarse-Graining of Particles .......................................................................................116 5.2 Conservative Interaction Parameters ...........................................................................118 Chapter 6 EMIm-Based Electrolytes ...................................................................................120 6.1 Structure of EMImAlX4 and EMImAl2X7 with MgX2 (X = Cl, I) ..............................120 6.2 Vibrational Spectra of the Electrolytes ........................................................................124 6.2.1 (EMImCl/(AlCl3)1.5)/(-MgCl2)x System .............................................................124 6.2.2 (EMImAlI4)/(-MgI2)x System .............................................................................127 Chapter 7 Mixed Acid Electrolytes ......................................................................................129 7.1 Hydration of the Acids ................................................................................................129 7.2 Proton Transfer ............................................................................................................144 Chapter 8 Conclusions ..........................................................................................................152 References ................................................................................................................................158 Appendices ...............................................................................................................................196 ix Appendix A Computer Programs and Routines ................................................................197 Appendix B List of Presentations and Publications ..........................................................216 Appendix C Reprints of Publications ................................................................................218 Vita ............................................................................................................................................358 x LIST OF TABLES Table 3.1. Average vanadium-oxygen distances for optimized structures shown in Figure 3.1.a ...............................................................................................................................................61 Table 3.2. Binding energies of water ligands to central vanadium cations calculated at B3LYP/6- 311G** and B3LYP/6-311+G** (in parenthesis) level of theory (all in kcal/mol). .....................63 Table 3.3. Calculated Gibbs energy of hydration of vanadium cations using the quasi-chemical approach compared with available experimental information (all in kcal/mol). ...........................64 Table 3.4. Average vanadium-oxygen distances for optimized structures shown in Figures 3.3 to 3.6.a ................................................................................................................................................70 Table 3.5. Computed atomic partial charges on vanadium, oxygen, and water molecules in the primary hydration shell using ChelpG scheme. .............................................................................77 Table 3.6. Selected bead types, their structures, and volumes. .....................................................84 Table 3.7. Calculated conservative interaction parameters. ..........................................................86 Table 3.8. Bonded interaction parameters obtained by fitting to the structures from classical MD simulations. ....................................................................................................................................87 Table 4.1. Selected bead types, their chemical structures, and volumes. ......................................96 Table 4.2. Calculated conservative interaction parameters. ..........................................................97 Table 4.3. Bonded interaction parameters for simulating the SEBS triblock copolymer and hydrated AEMs. The first four rows define the bonded parameters of the SEBS and all for the AEM. ..............................................................................................................................................99 Table 4.4. Diffusivities of W, OH−, and Cl− beads at various hydration levels for the hydrated AEMs at the alkaline and chloride forms. Diffusivities are in cm2/s × 104. ...............................114 xi Table 5.1. Selected bead types, their structures, and volumes. ...................................................117 Table 5.2. Calculated conservative interaction parameters. ........................................................119 Table 7.1. The number of water molecules required to effect dissociation of a proton in the gas phase and SMD solvated systems. The acids are ordered based on the acid strength (i.e., pKa values and gas phase acidities). ...............................................................................................................132 xii LIST OF FIGURES Figure 1.1. Schematic illustration of a VRFB during the discharge and charge processes (taken from reference 10). (a) Discharging: the V2+ cation oxidizes to V3+ and produces electrons that are transferred via an external circuit. The hydronium ion is transferred from the anolyte to the catholyte and the VO2 + cation reduces to VO2+. (b) Charging: the V3+ cation reduces to V2+ and the VO2+ cation oxidizes to VO2 +. The hydronium ion is transferred in a direction opposite to the discharge process. ............................................................................................................................3 Figure 1.2. Schematic representation of a PEMFC (taken from reference 62). Hydrogen and oxygen gases are inserted at the anode and cathode, respectively. Hydrogen gas produces protons and electrons. An external circuit transfers the electrons to the cathode and the protons permeate through the membrane electrolyte. At the cathode the oxygen and protons react and produce water. ...............................................................................................................................................7 Figure 1.3. The chemical structure of Nafion ionomer. The values of y and x determine the molecular weight and equivalent weight (molecular weight per ionic group), respectively. ..........9 Figure 1.4. Schematic representation of an alkaline fuel cell with a polymer electrolyte membrane. Hydrogen and oxygen gases insert at the anode and cathode, respectively. Hydrogen gas reacts with hydroxide at the anode generating water and electrons. The electrons are transferred through an external circuit to the cathode, where oxygen reacts with water to generate hydroxide ions. ..............................................................................................................................13 Figure 1.5. Chemical structures of common functional groups used in AEMs. (a) pyridinium; (b) quaternary ammonium; (c) phosphonium; (d) sulfonium; (e) guanidinium; and (f) imidazolium (taken from reference 98)...............................................................................................................14 xiii Figure 1.6. The chemical structures of: (a) the SEBS triblock copolymer; and (b) the synthesized AEM ionomer obtained by functionalization of the SEBS copolymer with quaternary ammonium groups (taken from reference 112). The two end blocks of the SEBS consist of polystyrene, a rigid polymer, while the midblock is hydrogenated polybutadiene, an elastomer. ................................15 Figure 2.1. Classification of several computational methods in chemistry and materials science as a function of length and time scales. The methods at larger length scales probe longer times but lack any small scale information....................................................................................................22 Figure 2.2. Coarse-graining of Nafion, a PFSA ionomer. Different mesoscopic particles are defined to represent chemically distinct subunits of the polymer (taken from reference 173)......23 Figure 2.3. Illustration of the conformational coordinate notation, 𝑅𝑛 (taken from reference 195). ...............................................................................................................................................47 Figure 3.1. First hydration shell structures of vanadium ions at the B3LYP6-311+G**/SMD level of theory: (a) V2+; (b) V3+; (c) VO2+; (d) VO2 +, octahedral; and (e) VO2 +, bipyramidal. Color scheme: oxygen (red), hydrogen (white), vanadium (purple). Two-tone solid lines present covalent bonds and dashed lines show coordinate bonds and hydrogen bonds (hydrogen bond criterion 1.1 Å < O…H < 2.5 Å). ........................................................................................................................61 Figure 3.2. Quasi-chemical contributions of the solvation free energy of aqueous V2+, V3+, VO2+, and VO2 + versus n number of water molecules inside the inner-shell structure. The calculations performed at the B3LYP/6-311+G(2d,p)//B3LYP/6-311G** level of theory and the outer-shell contributions are assessed using the SMD dielectric continuum model. Triangles represent inner- shell contributions and squares are outer-shell contributions to the total free energy of solvation. Circles are the total free energy of solvation predicting n = 6 as the most probable inner-shell occupancy for V2+ and V3+, n = 4 for VO2+ and n = 2 for VO2 +. ...................................................66 xiv Figure 3.3. Fully optimized configurations at B3LYP/6-311G**/SMD level of theory for the hydrated V2+ cation in proximity to: (a) H2SO4; (b) HSO4 −; (c) SO4 2−; (d) CF3SO3H; and (e) CF3SO3 −. The atomic color scheme is: oxygen (red), hydrogen (white), vanadium (purple), sulfur (yellow), carbon (grey), and fluorine (cyan). Distances are represented in Å, two-tone solid lines present covalent bonds, and dashed lines show coordinate bonds and hydrogen bonds (hydrogen bond criterion 1.1 Å < O…H < 2.5 Å). ...........................................................................................69 Figure 3.4. Fully optimized configurations at B3LYP/6-311G**/SMD level of theory for the hydrated V3+ cation in proximity to: (a) H2SO4; (b) HSO4 −; (c) SO4 2−; (d) CF3SO3H; and (e) CF3SO3 −. Color scheme as described in Figure 3.3. ......................................................................71 Figure 3.5. Fully optimized configurations at the B3LYP/6-311G**/SMD level of theory for the hydrated VO2+ cation in proximity to: (a) H2SO4; (b) HSO4 −; (c) SO4 2−; (d) CF3SO3H; and (e) CF3SO3 − . Color scheme as described in Figure 3.3. .....................................................................73 Figure 3.6. Fully optimized configurations at B3LYP/6-311G**/SMD level of theory for the hydrated VO2 + cation in proximity to: (a) H2SO4; (b) HSO4 −; (c) SO4 2−; (d) CF3SO3H; and (e) CF3SO3 − . Color scheme as described in Figure 3.3. .....................................................................75 Figure 3.7. The six degrees of freedom that contribute different configurations to the two objects (beads). The atomic structures of the beads 𝑖 and 𝑗 are arbitrary. The angles can vary as 0° ≤ 𝜃 ≤ 180° and 0° ≤ 𝜑 ≤ 360° for each configuration. The lower limit for the 𝑟𝑖𝑗 is when the beads start to overlap and the upper limit is the cut off distance 𝑟𝑐. ........................................................81 Figure 3.8. Distribution of: (a) the bond distance between connected A beads; (b) the angle between three consecutive A beads; and (c) the distance between two neighboring C beads in each Nafion chain. In order to compare atomistic chain structures from MD simulations to those of DPD, each Nafion chain was dissected into fragments equivalent to the selected beads and the center of xv mass of each fragment was considered as the position of the representing bead. The solid lines are the results from the MD simulation. The dashed lines are the corresponding DPD system with conservative interaction parameters from Table 3.7 and the default set of bonded interactions (𝑓𝑖 𝐵𝑜𝑛𝑑𝑒𝑑 = 4.0𝑟𝑏). The dash-dots are the corresponding DPD system with both bonded and non- bonded interactions from Table 3.7 and 3.8. .................................................................................88 Figure 3.9. RDFs of vanadium and C− beads for hydrated Nafion systems at different hydration levels and EW = 1115. (a) Vn+–Vn+ RDFs at λ = 6; (b) Vn+–Vn+ RDFs at λ = 18; (c) C−– C− RDFs at λ = 6; and (d) C−– C− RDFs at λ = 18. Color code: violet-V2+; green-V3+; blue-VO2+; yellow- VO2 +; black-hydrated Nafion. ........................................................................................................90 Figure 3.10. Radial distribution functions of the water beads at EW = 1115 and two different hydration levels: (a) 𝜆 = 6; and (b) 𝜆 = 18. The insets are magnification of the first peak in each graph. The solid line is for the system of hydrated Nafion, the dashed line has absorbed V2+, and the dotted line has absorbed V3+. ...................................................................................................91 Figure 3.11. Calculated vanadium and water diffusivities at different hydration levels in hydrated Nafion with an EW = 1115. (a) Vn+ beads; and (b) W beads. The V3+ and VO2 + beads show the lowest and the highest diffusivities, respectively. The effect of vanadium cations on the diffusivity of the water beads is insignificant. .................................................................................................92 Figure 3.12. RDFs and diffusivities for hydrated Nafion systems at different equivalent weights and λ = 12. (a) W–W RDFs; (b) VO2+– C− RDFs; (c) VO2+– VO2+ RDFs; and (d) diffusivities of W and VO2+ beads. The cation-anion interactions are increased by increasing the EWs, whereas the vanadium diffusivity is decreased. ...........................................................................................93 Figure 4.1. Coarse-grained model of: (a) SEBS triblock copolymer; and (b) hydrated AEM in alkaline form. Each circle represents the selected chemical group or collection of atoms for a xvi particular bead and is distinguished with a different color. The midblock of SEBS is modeled as polyethylene instead of the poly(ethylene-co-butylene) for simplicity (see Figure 1.5). ..............95 Figure 4.2. Simulated morphology of the SEBS triblock copolymer. Only the polystyrene end blocks, i.e., the Bs and Ph beads, are shown for clarity. (a) and (b) show two different views of the hexagonal cylindrical morphology. 𝑅𝑐𝑦𝑙𝑖𝑛𝑑𝑒𝑟 is the cylinder radius and 𝑑𝑐𝑦𝑙𝑖𝑛𝑑𝑒𝑟 is the intercylinder distance. The orange and mauve spheres illustrate the Bs and Ph beads, respectively. The simulation box is shown in blue. ..........................................................................................100 Figure 4.3. The polystyrene (i.e., BsPh – BsPh) RDF of the SEBS triblock copolymer. The position of the second point with the unit function value is twice the cylinder radius (2𝑅𝑐𝑦𝑙𝑖𝑛𝑑𝑒𝑟) and the position of the next peak is the intercylinder distance (𝑑𝑐𝑦𝑙𝑖𝑛𝑑𝑒𝑟). ...............................101 Figure 4.4. Simulated morphology of the hydrated AEM in the alkaline form at 𝜆 = 8. Each figure shows a specific bead type for clarity: (a) all beads; (b) Bs, Ph, TMM, and TMA+; (c) Bm; (d) TMA+; (e) OH−; and (f) W. (a) shows a strong phase separation between the functionalized polystyrene end blocks and the midblocks that are shown separately in (b) and (c). (d) and (e) illustrate the distribution of ionic species and (d) shows that the majority of water are in the polystyrene phase with a few distributed in the midblock phase. Color code: Bs and Bm (orange); Ph (mauve); TMM (green); TMA+ (purple); OH− (cyan); and W (blue). ....................................102 Figure 4.5. RDFs of the hydrated AEM in the alkaline form at 𝜆 = 8 . The inset shows a magnification at short distances. Color scheme: BsPh – BsPh (blue); BsPh – Bm (red); Bm – Bm (orange); and TMA+ – TMA+ (purple). .......................................................................................103 Figure 4.6. Simulated morphology of the hydrated AEM in the alkaline form for 𝜆 = : (a) 4; (b) 8; (c) 12; (d) 16; and (e) 20. All beads are shown and the color scheme is the same as in Figure 4.4. The ionomer phase separates into hydrophilic and hydrophobic phases at all hydration levels. xvii As the hydration increases the morphology changes from perforated and interconnected lamellae (𝜆 = 4 and 8) to perfect lamellae (𝜆 = 12) and then to disordered bicontinuous domains (𝜆 = 16 and 20). Domains exclusively consisting of water are formed within the hydrophilic phase at the two highest water contents only. ..................................................................................................105 Figure 4.7. Polystyrene RDFs for the hydrated AEM in the alkaline form at five levels of hydration. Increasing the water content increases the intensity of the first peak and decreases that of the second. The distance at which the RDF reaches one increases through hydration indicating that the hydrophilic phase swells. ...............................................................................................106 Figure 4.8. Cross sections of the morphologies (shown in Figures 4.6 (d) and (e)) in the xy, yz, and xz planes (top to bottom). (a) 𝜆 = 16; and (b) 𝜆 = 20. The exclusive water domains are within the hydrophilic phase and larger at the higher level of hydration. The cross sections reveal that the water domains are three dimensional. ..........................................................................................107 Figure 4.9. W – W RDFs for the hydrated AEM in the alkaline form at various degrees of hydration. Clearly, increasing the hydration intensifies the peaks. The distance at which the RDF reaches one is the same for 𝜆 = 8 and 12 and increases for 𝜆 = 16 and 20. This is due to the formation of domains that exclusively contain water. There are not sufficient W beads in the simulation box at 𝜆 = 4 to form a smooth curve and therefore this data is not shown. ..............108 Figure 4.10. (a) TMA+ – OH− RDFs for the hydrated AEM in the alkaline form at various levels of hydration. (b) An expanded view of the first peak seen in (a). (c) The integrated number of the OH− beads around the TMA+ beads. There is one OH− in the immediate vicinity of each TMA+. At a higher degree of hydration there are fewer numbers of OH− around TMA+ at r > 5 Å. ......109 Figure 4.11. Simulated morphology of the hydrated AEM in the chloride form at 𝜆 = 8. Each figure shows a specific bead type for clarity: (a) all beads; (b) TMA+; (c) Cl−; and (d) W. The AEM xviii ionomer shows phase separation similar to the alkaline form. Color code: Bs and Bm (orange); Ph (mauve); TMM (green); TMA+ (purple); Cl− (shamrock green); and W (blue). .........................110 Figure 4.12. Simulated morphology of the hydrated AEM in the chloride form for 𝜆 = : (a) 4; (b) 8; (c) 12; (d) 16; and (e) 20. All beads are shown and the color scheme is the same as in Figure 4.11. The ionomer phase separates into hydrophilic and hydrophobic phases at all hydration levels. Similar to the alkaline form, the morphology changes by increasing the water content and exclusive domains of water are formed within the hydrophilic phase at the two highest values of 𝜆. ........111 Figure 4.13. W – W RDFs for the hydrated AEM in the alkaline (solid lines) and chloride (dashed lines) forms at various levels of hydration. Changing the anion from OH− to Cl− results in negligible changes at 𝜆 = 8 and 12, however it significantly increases the peaks at 𝜆 = 16 and 20. The distance at which the RDF reaches one substantially increases at 𝜆 = 16. ..........................112 Figure 4.14. Cross sections of the hydrated AEM morphologies in the: xy (top), yz (middle), and xz (bottom) planes at 𝜆 = 16. (a) Alkaline form; and (b) chloride form. Changing the anion from OH− to Cl− results in larger exclusive domains of water. ............................................................113 Figure 5.1. The repeat unit of the SDAPP ionomer. The number of sulfonate groups per monomer depends on the level of sulfonation which is controlled by the reaction time (taken from reference 74). ...............................................................................................................................................116 Figure 5.2. The coarse-graining scheme for the SDAPP ionomer. (a) The monomer without sulfonation; and (b) two sulfonate groups per monomer. Different ratios of (a) and (b) may be combined to represent sulfonation levels ranging from 0 to 2. Different bead types are shown by colored ovals. Color code: P1 (red), P2 (blue), P3 (green). .........................................................118 Figure 6.1. Optimized geometries of the ionic liquid complexes determined at the B3LYP/6- 311G** level of theory. (a) EMImAlCl4; (b) EMImAlCl4 with MgCl2; (c) EMImAlCl4 with xix (MgCl2)2; (d) EMImAl2Cl7; (e) EMImAl2Cl7 with MgCl2; and (f) EMImAl2Cl7 with (MgCl2)2. Color code: chlorine (green), aluminum (pink), magnesium (yellow), nitrogen (blue), hydrogen (white), carbon (grey). Distances between chlorine and hydrogen are only shown when less than 3 Å. ...............................................................................................................................................122 Figure 6.2. Optimized geometries of the ionic liquid complexes determined at the B3LYP/6- 311G** level of theory. (a) EMImAlI4; (b) EMImAlI4 with MgI2; (c) EMImAlI4 with (MgI2)2; (d) EMImAl2I7; (e) EMImAl2I7 with MgI2; and (f) EMImAl2I7 with (MgI2)2. Color code: iodine (orange), aluminum (pink), magnesium (yellow), nitrogen (blue), hydrogen (white), carbon (grey)............................................................................................................................................123 Figure 6.3. Electrostatic potential map on the total electron density isosurface. (a) [AlCl4MgCl2] −; (b) [Al2Cl7MgCl2] −; (c) [AlCl4(MgCl2)2] −; and (d) [Al2Cl7(MgCl2)2] − complexes. These maps show a more negative electrostatic potential and therefore more electron rich area on MgCl2 rather than aluminium chloride. Color code: chlorine (green), aluminum (pink), magnesium (yellow). ...................................................................................................................124 Figure 6.4. Vibrational spectra of ([EMIm]Cl/(AlCl3)1.5)/(-MgCl2)x electrolytes with 0 ≤ x ≤ 0.20. (a) Raman spectra in the region between 100 and 480 cm−1, the intensity of the band corresponding to the stretching of the monomeric species AlCl4 − (blue) increases at higher concentrations of the Mg salt, whereas the intensity of the band attributed to the stretching of the dimeric species Al2Cl7 − (yellow) decreases; and (b) far IR spectra, the intensities of the bands corresponding to the vibrations of polymeric MgCl2 species (green) and Mg–Cl–Al stretching modes (blue and yellow) increase at higher concentrations of the Mg salt. .........................................................................126 Figure 6.5. Vibrational spectra of (EMImAlI4)/(δ-MgI2)x electrolyte with 0 ≤ x ≤ 0.023. (a) Raman spectra with band assignments based on DFT calculations, the intensity of the band assigned to the xx stretching of AlI4 − (blue) increases with the δ-MgI2 concentration while that of Al2I7 − (orange) decreases; and (b) far IR spectra. .................................................................................................127 Figure 7.1. Fully optimized structures of CF3SO3H with from zero to five water molecules at the B3LYP/6-311G** level of theory: (a) gas phase structures, proton dissociation occurs after three water molecules are added; (b) structures with the SMD model, proton dissociates after only one water molecule is added. The atomic color scheme in Figures 7.1 to 7.10 is: oxygen (red), hydrogen (white), sulfur (yellow), carbon (gray), fluorine (cyan), nitrogen (blue), phosphor (orange). ....131 Figure 7.2. Fully optimized structures of FSO3H with from zero to five water molecules at the B3LYP/6-311G** level of theory: (a) gas phase structures, proton dissociation occurs after three water molecules are added; (b) structures with the SMD model, proton dissociates after only one water molecule is added. ..............................................................................................................133 Figure 7.3. Fully optimized structures of CH3SO3H with from zero to five water molecules at the B3LYP/6-311G** level of theory: (a) gas phase structures, proton dissociation occurs after three water molecules are added; (b) structures with the SMD model, proton dissociates after only two water molecules are added. ..........................................................................................................134 Figure 7.4. Fully optimized structures of H2SO4 with from zero to five water molecules at the B3LYP/6-311G** level of theory: (a) gas phase structures, proton dissociation occurs after four water molecules are added; (b) structures with the SMD model, proton dissociates after three water molecules are added. In all structures, the water molecules bind to one side of the acid. The only exception is the gas phase with three water molecules. ...............................................................136 Figure 7.5. Fully optimized structures of HNO3 with from zero to five water molecules at the B3LYP/6-311G** level of theory: (a) gas phase structures, proton dissociation occurs after five xxi water molecules are added; (b) structures with the SMD model, proton dissociates after only three water molecules are added. ..........................................................................................................137 Figure 7.6. Fully optimized structures of HF with from zero to five water molecules at the B3LYP/6-311G** level of theory: (a) gas phase structures; (b) structures with the SMD model. No proton dissociation occurs in agreement with the weak acidity of HF. .................................138 Figure 7.7. Fully optimized structures of H3PO2 with from zero to five water molecules at the B3LYP/6-311G** level of theory: (a) gas phase structures; (b) structures with the SMD model. No proton dissociation occurs, except in the gas phase with five added water molecules. ......140 Figure 7.8. Fully optimized structures of H3PO3 with from zero to five water molecules at the B3LYP/6-311G** level of theory: (a) gas phase structures; (b) structures with the SMD model. No proton dissociation occurs. .....................................................................................................141 Figure 7.9. Fully optimized structures of H3PO4 with from zero to five water molecules at the B3LYP/6-311G** level of theory: (a) gas phase structures; (b) structures with the SMD model. No proton dissociation occurs. .....................................................................................................142 Figure 7.10. Fully optimized structures of CF3SO2NHSO2CF3 with from zero to five water molecules at the B3LYP/6-311G** level of theory: (a) gas phase structures, proton dissociation occurs after three water molecules are added (nitrogen is not involved in the hydrogen bond network); (b) structures with the SMD model, proton dissociates after only two water molecules are added (nitrogen is involved in the hydrogen bond network). ................................................143 Figure 7.11. The relative energy profiles of PES scans for: CF3SO3H, FSO3H, and CH3SO3H. The gas phase results are shown in the left panel and the results with the SMD model are displayed in the right panel. The extremum plots in the gas phase correspond to the back transfer events indicating that the molecular acid is more stable than the ionized form. The gas phase results are xxii similar, because the proton dissociation occurs after three water molecules are added. The SMD results for CF3SO3H and FSO3H are similar to each other and different to those for CH3SO3H, as the formers require only one water molecule to effect proton dissociation. The right panel scans do not show a back transfer of a proton, because the formed hydronium ion is stabilized by the continuum solvation model. .........................................................................................................146 Figure 7.12. The relative energy profiles of PES scans for: H2SO4, HNO3, and HF. The gas phase results are shown in the left panel and the results with the SMD model are displayed in the right panel. The gas phase results exhibit back transfer of a proton indicating that the molecular acid is more stable than the ionized form. The energy barrier to move a proton away from the acid is inversely related to the number of water molecules and the acid strength. The weakest acid, HF, has the highest energy barriers. The relative energy axis is plotted from 0 to 18 kcal/mol for all acids except for HF (from 0 to 35 kcal/mol)................................................................................148 Figure 7.13. The relative energy profiles of PES scans for: H3PO2, H3PO3, and H3PO4. The gas phase results are shown in the left panel and the results with the SMD model are displayed in the right panel. Only a back transfer or an increasing scan is observed, because the proton dissociation is not favored by any of these weak acids. The proton transfer energetics are similar for these acids in agreement with their similar acidity and pKa values listed in Table 7.1. .................................149 Figure 7.14. The relative energy profiles of PES scans for CF3SO2NHSO2CF3. The gas phase results are shown in the left panel and the results with the SMD model are displayed in the right panel. In the gas phase with one water molecule, the secondary proton is transferred to an oxygen atom of the acid. The gas phase structure of this acid with three water molecules does not involve hydrogen bonding with nitrogen as shown in Figure 7.10 (a). Hence the proton was transferred to an oxygen atom of the acid, i.e., the orange plot in the left panel shows OH distance rather than xxiii NH. The proton dissociation occurs after two water molecules are added to the SMD solvated structures, however the proton transfer energetics is negligible with one water molecule. ........150 1 Chapter 1 Introduction Recent interest in replacing fossil fuels with renewable green energy sources such as solar and wind has significantly increased as a result of the perceived consequences to the environment. One of the limiting factors for the utilization of renewable sources in either the electric distribution grid or smart grid is the intermittent nature of these sources which requires the employment of electrical energy storage devices.1 The conversion of energy by clean, efficient, and sustainable means such as fuel cells is another challenge in the replacement of fossil fuels. Lead acid, Ni-Cd, lithium- ion/air, and magnesium batteries along with fuel cells and redox flow batteries constitute the majority of the current energy storage and conversion technologies.2 The specific systems of interest in this body of work are the all-vanadium redox flow battery, the alkaline fuel cell, and the secondary magnesium battery. 1.1 Redox Flow Batteries Redox Flow Batteries (RFBs) are among the promising large-scale (grid size) electrical energy storage systems and have received attention due to their superior properties. These include: long lifespan; simple installation; low maintenance cost; the possibility of instant charging through the replacing of the electrolytes; and the ability to repeatedly store and convert electrical energy to chemical energy and vice versa.1 A RFB stores energy in two electrolyte solution tanks containing different redox couples. The liquid electrolyte solutions flow through the electrochemical cell, where electrochemical reactions occur. The two electrolytes are separated by an ion exchange membrane within the cell. The energy storage capacity is determined by the size of the electrolyte tanks, and the power is a dependent upon the size of the cell stacks. Thus, the power and energy of a RFB are decoupled.3 Various 2 redox couples have been employed in RFBs, including but not limited to: Cr2+/Cr3+ vs. Fe2+/Fe3+;4- 5 S2 2−/S4 2− vs. Br−/Br2; 6-7 Zn/Zn2+ vs. Br−/Br2; 8 and V2+/V3+ vs. V4+/V5+.9 The last mentioned in this list is referred to as the vanadium redox flow battery as it uses all-vanadium electrolytes. 1.1.1 Vanadium Redox Flow Batteries (VRFBs) The VRFB demonstrates the greatest commercial prospects among the many reported RFBs, due to its high energy efficiency, long cycle-life, and reversibility.3,9 Having all the active species from one element at different oxidation states, makes the regeneration of the electrolytes relatively easy. The VRFB consists of a V2+/V3+ sulfate solution as the negative electrolyte (i.e., anolyte) and a VO2+/VO2 + sulfate solution as the positive electrolyte (i.e., catholyte) and are separated by an ion selective membrane that permits the transport of protons. A schematic of the VRFB system is illustrated in Figure 1.1. Figure 1.1 (a) shows the VRFB during the discharge process where the VO2 + is reduced to VO2+ and the V2+ cation is oxidized to V3+. The hydronium ion is transferred from the anolyte to the catholyte to complete the circuit and generate electricity. The same reactions take place during the charging (see Figure 1.1 (b)) but in the opposite direction, i.e., the VO2 + and V2+ cations are oxidized and reduced, respectively and the hydronium ion is transferred from the catholyte to the anolyte. The half-cell reaction during charging for the negative electrolyte is: 𝑉3+ + 𝑒−⟶ 𝑉2+ ; (1.1) and for the positive electrolyte is: 𝑉𝑂2+ +𝐻2𝑂 ⟶ 𝑉𝑂2 + + 2𝐻+ + 𝑒− . (1.2) Hence, the total cell reaction may be written as: 𝑉𝑂2+ +𝐻2𝑂 + 𝑉 3+⟶ 𝑉𝑂2 + + 2𝐻+ + 𝑉2+ . (1.3) 3 Figure 1.1. Schematic illustration of a VRFB during the discharge and charge processes (taken from reference 10). (a) Discharging: V2+ is oxidized to V3+ and produces electrons that are transferred via an external circuit. The hydronium ion is transferred from the anolyte to the catholyte and the VO2 + cation is reduced to VO2+. (b) Charging: the V3+ cation is reduced to V2+ and the VO2+ cation is oxidized to VO2 +. The hydronium ion is transferred in a direction opposite to the discharge process. (a) (b) 4 The discharge process obeys the same reactions, but in the opposite direction. The wide difference in the reduction potentials of the two half-cells produces a standard voltage of 1.25 V, which is comparable to other types of RFBs.10 1.1.2 Electrolytes for VRFBs The energy density of the battery depends on the concentration of the electrolyte solutions. VRFBs have low energy density and therefore need large electrolyte tanks which limit their application to stationary storage.1,11 Hence, increasing the vanadium concentration results in an increase in the energy density of the battery. However, the concentration of the vanadium ions is limited by the stability of the electrolytes and precipitation of the four vanadium cations. The stability of the electrolytes is determined by the concentration of sulfuric acid, solution temperature, and the state of charge of the battery.12-14 The stability of VO2 + (in H2SO4) is compromised by V2O5 precipitation which occurs at high temperature (> 40 °C) and high vanadium concentration (> 2 M).15 Although, the VO2 + concentration can be increased by the addition of sulfuric acid to the solution, the increase in acid concentration is limited by the stability of the VO2+ species, as both VO2 + and VO2+ species are present in the positive electrolyte.14 In contrast to VO2 +, the V2+, V3+, and VO2+ solutions are stabilized by increasing the temperature and decreasing the sulfuric acid concentration,12,14 which makes the design of VRFB electrolytes challenging. Inorganic and organic agents have been added in an effort to stabilize the electrolytes.13,16-19 The stability of the vanadyl solution V(IV) at low temperatures can be improved through the addition of stabilizing agents, however the effect of the stabilizers must be considered on both the V(IV) and V(V) solutions as they co-exist in the positive electrolyte. Furthermore, since the crossover of stabilizing agents through the membrane is possible, their effects need to be considered on all four vanadium cations simultaneously.13,18 5 Under the high ionic strength condition of a VRFB, understanding the hydration of the vanadium cations and their interaction with the solvent and membrane is crucial in optimizing the design and modification of the electrolytes and membranes. Vanadium complexes have been subject to both experimental and theoretical investigations.20-28 However, the structure of hydrated vanadium ions, the oxidation state of the ions, and their interaction with sulfuric acid and/or a hydrated perfluorosulfonic acid (PFSA) ionomer are not understood.29 Hence, a molecular-level investigation on the interaction of vanadium cations in a VRFB electrolyte is warranted. Membranes used in VRFBs are classified as either cation exchange or anion exchange depending on the type of ion that is permeated. Both types of membranes must exhibit low permeability to the vanadium ions, high proton or sulfate conductivity, good thermal and chemical stability, substantial mechanical integrity, and low cost.3,11 Typically, the imperfect ion selectivity of the cation exchange membrane results in vanadium cation crossover that causes electrolyte contamination and capacity loss of the battery as well as a decrease in conductivity of the membrane.3,30-32 The diffusion of sulfuric acid into Nafion has also been observed.33-35 Hence, there is ongoing research on the development of new membranes and membrane treatments to decrease solvent and vanadium ion transport across the membrane.36-39 One of the most widely used cation exchange membranes for VRFBs is the PFSA membrane, Nafion.40 Nafion phase separates when hydrated into nano-dimensioned regions containing water associated with the sulfonic acid groups surrounded by the polytetrafluoroethylene backbone of the ionomer.41-45 Protons and vanadium cations permeate through the membrane (via the water domains). Although, the mobility of the vanadium cations is much lower than that of a proton, over time the crossover of vanadium cations results in significant capacity loss and electrolyte contamination.30,46-47 6 Understanding the diffusion of vanadium cations through the membrane and the morphology evolution, from a meso-scale perspective, is another important aspect in the design of a VRFB, due to impact on electrolyte contamination, battery capacity, and the hydration and conductivity of the membrane. The dynamics and morphological properties of hydrated Nafion have been simulated extensively,48-53 however the effect of the absorbed vanadium cations are not known. There is a recent molecular dynamics (MD) study on the effect of the absorbed V2+ and V3+ species on the hydrated Nafion,54 and more investigations are required to understand their effects. 1.2 Fuel Cells Fuel cells are an important sustainable energy conversion devices. These devices convert chemical energy into electrical energy through electrochemical reactions of a gaseous fuel (e.g. hydrogen or methanol) and an oxidant gas (e.g. oxygen from the air) producing water.55 Energy conversion via fuel cells is pollution free and more efficient than other conventional thermo-mechanical technologies.55-56 Fuel cells may be classified into five major groups based on the type of the electrolyte: alkaline fuel cells (AFCs), proton-exchange membrane fuel cells (PEMFCs), solid oxide fuel cells (SOFCs), phosphoric acid fuel cells (PACFs) and molten carbonate fuel cells (MCFCs).57 PEMFCs and AFCs have received significant attention due to their application and commercialization potential58-60 and are further described below. 1.2.1 Proton Exchange Membrane Fuel Cells (PEMFCs) The first PEMFC was designed and developed in 1955 at General Electric (GE).61 PEMFCs use a solid polymer electrolyte that results in safer and easier cell operation and handling. PEMFCs may be utilized in automotive vehicles due to their quick start up and high power density. Despite many advantages they are restricted to operating temperatures below 100 °C and also require expensive platinum catalysts.56 7 Figure 1.2. Schematic representation of a PEMFC (taken from reference 62). Hydrogen and oxygen gases are inserted at the anode and cathode, respectively. Hydrogen gas is oxidized producing protons and electrons. An external circuit transfers the electrons to the cathode and the protons permeate through the membrane electrolyte. At the cathode oxygen is reduced producing only water. A schematic design of the PEMFC is illustrated in Figure 1.2. Hydrogen and oxygen gases are introduced at the anode and cathode, respectively. Hydrogen gas produces protons and electrons. An external circuit transfers the electrons to the cathode and the protons permeate through the electrolyte (membrane). At the cathode the oxygen and protons react and produce water. The electrochemical reaction at the anode is called the hydrogen oxidation reaction (HOR) and is written as: 𝐻2⟶ 2𝐻+ + 2𝑒− ; (1.4) and at the cathode the oxygen reduction reaction (ORR) takes place: 2𝐻+ + 1 2 𝑂2 + 2𝑒 −⟶𝐻2𝑂 . (1.5) The overall cell reaction is: (Proton Exchange Membrane) 8 𝐻2 + 1 2 𝑂2⟶𝐻2𝑂 . (1.6) The ORR at the cathode is kinetically slow and requires relatively high platinum loading that results in high cell costs.58,60 1.2.2 Electrolytes for PEMFCs The solid polymer electrolyte is the central component of the PEMFC and is required to possess high proton conductivity, adequate mechanical strength, high chemical stability, low fuel gas permeability, and low cost.59 The membrane materials may be classified as follows:  Perfluorinated ionomers  Non-fluorinated ionomers  Acid-base complexes  Supported composites The following provides a brief description and several examples for each group. Perfluorinated Ionomers. The perfluorinated polymers functionalized by sulfonic acid groups, i.e., PFSA membranes, have been extensively studied and utilized in PEMFCs due to their high proton conductivity and substantial thermal and mechanical stability.59 The high proton conductivity of PFSAs is due to the hydration and consequent proton dissociation of the strong acid groups. Nafion, originally developed by DuPont, is the most widely used PFSA ionomer as the cation exchange membrane for both PEMFCs and VRFBs.40,56,63 Figure 1.3 shows the chemical structure of Nafion with the equivalent weight (EW) being the molecular weight per ionic group of the ionomer. The polytetrafluoroethylene (PTFE) backbone provides outstanding chemical stability in both basic and acidic environments. 9 Figure 1.3. The chemical structure of Nafion ionomer. The values of y and x determine the molecular weight and equivalent weight (molecular weight per ionic group), respectively. Despite many advantages, PFSA membranes have high costs, limited operating temperatures, permeability to methanol, complex water management, and environmental issues related to the recycling.64 These challenges have led to significant research into finding alternative materials either by modifying the PFSAs or synthesizing and developing novel polymeric and composite materials. One of the approaches to modifying the present PFSAs is to replace the sulfonic acid group with other protogenic groups. The bis[(trifluoromethane)sulfonyl]imide, CF3SO2NHSO2CF3, is one of the alternatives that shows higher proton conductivity at intermediate hydration levels when substituted for the sulfonic acid group in Nafion.65 This sulfonyl imide is one of the strongest superacids66 and may be used under more extreme conditions.67-68 Phosphonic acid group is another alternative owing to its good thermal and oxidative stability.69 Although phospho-acids (phosphoric, phosphonic, and phosphinic) are not as strong as sulfonic acids under certain conditions they may exhibit even higher proton conductivity. Phosphoric acid is an outstanding proton conductor as a neat liquid.70-72 If transferrable to phosphonic acid this may lead to improved conductivities over the PFSAs at very low hydrations and anhydrous conditions. Burton et al. prepared perfluorophosphonic acid ionomers that showed good thermal stability and equal or lower conductivities to that of Nafion control.73 10 Non-Fluorinated Ionomers. The non-fluorinated hydrocarbon polymers have received tremendous attention as low-cost alternatives to PFSAs.62,64 These membranes are typically made from aromatic or heterocyclic repeat units to enhance the stability at elevated temperatures. Polysulfones (PSU), poly(ether sulfone) (PES), poly(ether ketone)s (PEK) with varying number of ether and ketone functionalities (e.g., PEEK, PEKK, PEEKK), and polybenzimidazole (PBI) are some examples of this group. These polymers are inherently insulating and therefore need to be modified for application in PEMFCs. Several approaches have been developed to add the ion conducting functional groups to the polymers including: acid or base doping; direct sulfonation of the backbone; grafting of a sulfonated or phosphonated functional group; and synthesis from the sulfonated monomer.64 The characteristics of the final membrane depend on the chemical nature of the polymer backbone, the molecular weight of the polymer, the distribution of the molecular weight, and the degree of functionalization. A competitive hydrocarbon-based ionomer is a sulfonated Diels-Alder poly(phenylene) (SDAPP) polymer.74 SDAPP can achieve comparable conductivity to Nafion and has been shown to outperform it when used in VRFBs.75-77 Transmission electron microscopy (TEM) imaging shows that the ionic domains in SDAPP are roughly 0.5 nm in diameter compared to Nafion where the average dimension of the water domains are approximately 2 nm when fully hydrated.76 It is important to understand the factors controlling the size of the ionic domains in the hydrated ionomer in order to enhance both ionic conductivity and selectivity. Acid-Base Complexes. This class of membranes incorporates acid components into an alkaline polymer base to promote proton conduction without humidification at high temperatures.59,62 The phosphoric acid-doped polybenzimidazole (PBI/H3PO4) is a successful example showing high proton conductivity and excellent mechanical properties at temperatures up to 200 °C.78 Acid-base 11 polymer blends are also promising materials exhibiting low water uptake, high proton conductivity, good mechanical and thermal stability. The improved properties of these blends are related to the acid-base interactions between the polymer chains, specifically, ionic interactions and hydrogen bonding.79 The polymer blend of sulfonated poly(ether ether ketone) as the acidic compound and polybenzimidazole as the basic compound, i.e., sPEEK/PBI, is an example with good proton conductivity and excellent thermal stabilities.80 Supported Composites. These membranes are obtained by combining polymers with inorganic materials; either the polymer or the inorganic component or both may be ion conducting. Composites incorporate the advantages of both polymers (e.g., flexibility and ionic conductivity) and inorganic materials (e.g., thermal and mechanical stability) plus the synergistic effects due to their interactions.81 The composite membranes offer higher thermal stability, enhanced proton conductivity, water retention at high temperatures, and mechanical support when compared to the pure polymer membranes.62 Metal oxides, metal phosphates, phosphosilicates, porous inorganic materials (e.g., zeolite, activated carbon, and metal organic framework (MOF)), and carbon nanotubes are some of the common inorganic materials used. The preparation technique controls the final properties of the composites because it directly affects the particle dispersion (determines the isotropy) and interfacial adhesion (relates to the synergism).81 The following preparation methods: infiltration; recasting; and electrospinning may be utilized depending on the compositions of the mixing components and the desired characteristics. The first two methods incorporate inorganic particles into the polymer by swelling or dissolution in a solvent. The last technique involves electrospinning of a non-conductive or less conductive phase into a mat and then adding a highly proton conducting component.82-85 The 12 inverse, i.e., electrospinning a conductive phase and then enforcing it with a non-conductive material has also been done.86-89 1.2.3 Alkaline Fuel Cells (AFCs) The AFC is one of the oldest fuel cells and was first described by Reid in 1903.90-91 It uses an aqueous electrolyte solution of potassium hydroxide (KOH), because it is the most conducting of all alkaline hydroxides.57 Hydrogen gas reacts with hydroxide anions at the anode generating water and electrons. The electrons are transferred through an external circuit to the cathode, where oxygen reacts with water to generate hydroxide ions.57 The half-cell reaction at the anode may be written as: 𝐻2 + 2𝑂𝐻 −⟶ 2𝐻2𝑂 + 2𝑒 − ; (1.7) and at the cathode: 𝑂2 + 2𝐻2𝑂 + 4𝑒 −⟶ 4𝑂𝐻− . (1.8) The overall cell reaction is the same as in a PEMFC, i.e.: 2𝐻2 + 𝑂2⟶ 2𝐻2𝑂 . (1.9) Some of the advantages of the AFC over other fuel cells include: lower operating temperature (≈ 23 – 70 °C), higher reaction kinetics at the electrodes resulting in higher cell voltages, and much lower noble metal catalyst loading. The ORR in AFC is more facile than in PEMFC due to the basic environment. Thus, AFC may operate with little or no platinum catalyst reducing the cost of the cell. It has been shown that non-precious metals (silver, nickel, and cobalt)92-93 and transition metal oxides94-95 may successfully be used as the catalyst. One of the main drawbacks in the AFC is the sensitivity of the KOH solution to CO2. 96 This necessitates the use of a pure feed oxidant stream instead of air. It is difficult to purify the fuel stream and consequently limits the commercialization of AFCs. The KOH solution may be 13 replaced with an anion-conducting polymer electrolyte.97 This relatively new design has an anion exchange membrane (AEM) functioning as both the separator and the conductive support between the electrodes and is shown in Figure 1.4. No solid crystals of the metal carbonate are formed that would block the gas diffusion electrodes, because there are no mobile cations (K+) and hence, the efficiency and lifetime of AFC are improved. Figure 1.4. Schematic representation of an alkaline fuel cell with a polymer electrolyte membrane. Hydrogen and oxygen gases are introduced at the anode and cathode, respectively. Hydrogen gas reacts with hydroxide at the anode generating water and electrons. The electrons are transferred through an external circuit to the cathode, where oxygen reacts with water to generate hydroxide ions. 1.2.4 Electrolytes for Anion Exchange Membrane Fuel Cells AEMs must have high ionic conductivity (but not electron conducting), low gas permittivity, and high thermal and mechanical stability.57,98 They are generally polymers functionalized by cationic groups capable of conducting hydroxide ions (OH−). A few anion-exchange groups are depicted in Figure 1.5. Benzyl trimethyl quaternary ammonium (Figure 1.5 (b) with R = –CH3) is one of the most widely used cation exchange groups in AEMs.99 14 Figure 1.5. Chemical structures of common functional groups used in AEMs. (a) pyridinium; (b) quaternary ammonium; (c) phosphonium; (d) sulfonium; (e) guanidinium; and (f) imidazolium (taken from reference 98). The conductivity of AEMs is generally lower than PEMs, because: (1) the hydroxide ion has a lower mobility than that of the proton;100 and (2) some AEMs lack the microphase separation between the hydrophilic and hydrophobic phases that facilitate ion transport.60,101-102 A main problem with AEMs is the limited chemical and thermal stability under highly corrosive alkaline conditions. The cationic groups may be decomposed by OH− that results in degradation of the polymer, reduced ion exchange capacity, and lower conductivity of the membrane. It has been determined that the stability of an AEM strongly depends on the chemical nature of the cationic functional groups and the main chain structure.60 The alkaline stability may be improved by using stable cation structures,103-107 imposing steric hindrance and electronic effects,106,108-109 and crosslinking.110-111 (a) (b) (c) (d) (e) (f) 15 Several chemically stable AEMs using the triblock copolymer: polystyrene-b-poly(ethylene-co- butylene)-b-polystyrene (SEBS), functionalized by various benzyl- and alkyl-substituted quaternary ammonium groups have been recently developed.112-114 The chemical structures of the SEBS and SEBS functionalized by a quaternary ammonium group (SEBS-QA) are shown in Figure 1.6. The SEBS backbone provides good chemical and mechanical stability to the membrane due to its robust structure, high molecular weight, and elastomeric units. It is known that the morphology of the SEBS copolymer depends on the block compositions.115-118 Hence, the phase- separated morphology of SEBS-based AEMs may be tailored to achieve high anion conductivity by varying the block composition in the SEBS copolymer. Figure 1.6. The chemical structures of: (a) the SEBS triblock copolymer; and (b) the synthesized AEM ionomer obtained by functionalization of the SEBS copolymer with quaternary ammonium groups (taken from reference 112). The two end blocks of the SEBS consist of polystyrene, a rigid polymer, while the midblock is hydrogenated polybutadiene, an elastomer. Fundamental understanding of transport in AEMs is still in its infancy,119 in contrast to PEMs which have been systematically studied for many years.120 Quantum chemical calculations121-123 and molecular dynamics simulations124-127 have previously been used to determine the structure and stability of the anion exchange functional groups as well as the transport mechanism of hydroxide in AEMs, similar to what has been done for PEMs.128-132 These techniques provide (a) (b) 16 valuable information at the atomic and molecular scales, however they are inadequate to explore the morphological aspects due to the length and time scale limitations. Meso-scale simulations may provide a relative understanding of the morphological features of AEMs as they previously did for PEMs.48-49,133-134 The effect of other anions including chlorine and bromine may also be studied, because the membranes are usually prepared and characterized in a halogen precursor form.112,114 1.3 Secondary Magnesium Batteries Secondary batteries are of great interest to industry as they may offer solutions to the current technical and cost barriers as well as the safety and environmental concerns associated with portable energy storage.2 Magnesium based batteries are deemed particularly promising if advancements in materials and manufacturing technologies are achieved.2,135-136 These batteries show advantages over other competing types (i.e., lithium and sodium ion batteries) owing to the abundance of Mg in the earth’s crust and the issues concerning disposal and waste. They also have a broad operating temperature range and their capacity loss over time is quite small. The magnesium cation is, of course, divalent thereby providing twice the charge per ion as monovalent ions such as lithium. This results in a theoretical volumetric capacity of 3832 mAh/cm−3 (in contrast to Li which has a capacity of 2061 mAh/cm−3).135 It is also more stable and less chemically reactive than either Li or Na, and consequently, large scale magnesium electrical energy storage devices are simpler to design.2 A practical magnesium rechargeable battery was pioneered by Aurbach et al. with an energy density of 60 Wh/kg which was operated for more than 2000 cycles with little fade.137 Proper design and architecture may lead to even higher gravimetric energy densities of 150 − 200 Wh/kg and operating voltages between 2 and 3 volts.2,138 17 The anode is usually made of commercially pure magnesium (99.9% or higher) in different forms such as foil, ribbon, or coupon.2,137 At the cathode, where the electrons are transferred to the host material, Mg2+ ions must be inserted and extracted from a solid host to the electrolyte solution without causing any electrode destruction. Many materials fail to serve as the cathode for magnesium batteries due to the difficulty of intercalating magnesium ions.137 Cathodes made of chevrel phases (CPs), MxMo6T8 (M = metal, T = S, Se, Te) a unique class of host materials, exhibit fast and reversible intercalation.137,139 Other cathodes have been developed based on V2O5, MoO3, transition metal oxide/sulfide/boride, and NASICON-type materials.2 Recently, Muldoon and coworkers developed a new class of cathodes based on the conversion cathode reaction of metallic magnesium with elemental sulfur, referred to as the magnesium/sulfur conversion cathode.136,140 The combination of a magnesium anode and a sulfur cathode can theoretically deliver a high energy density (~ 4000 Wh/l). However, the greater challenge is to identify suitable electrolytes which are compatible with sulfur. A critical roadblock towards practical Mg-based energy storage is the lack of reversible electrolytes that are safe and electrochemically stable. The choice of the electrolyte and its properties also affects the class of cathodes utilized. Magnesium has a unique electrochemistry which prohibits its reversible deposition in aprotic solvents due to the formation of a passivation layer which prohibits conduction of Mg2+ ions.141 Ethereal solutions of Grignard reagents and magnesium organohaloaluminates may surpass this problem.139,142-143 However, these solutions are hazardous due to their volatility and flammability. Existing electrolytes136,144-146 based on liquid solvents are inadequate for meeting the needs of functional devices for portable electronics and transportation applications due to their high volatility and corrosiveness.135,138 The volatility issue may be relieved by using polymers, however, 18 the possibility of synthesizing high-performance polymer electrolytes with the  and  forms of MgCl2 had been ruled out due to the high lattice energy of the salt.147 MgCl2 may be prepared in a form called -MgCl2 that is characterized by high crystallographic disorder, reactivity, and solubility.148 The non-conventional properties of this salt are due to the presence of a metastable nanoribbon or polymeric structure with concatenating MgCl2 repeating units, in which the Mg atoms are bridged via chloride bridges. The preparation of polymer electrolytes using this salt has yielded Mg2+-conducting materials with conductivities as high as 10-4 S·cm-1 at room temperature. These studies lead to the first Mg-polymer batteries.149 The use of the -form of MgCl2 with ionic liquids as the electrolyte is a significant milestone in the search for high-performance stable electrolytes. 1.3.1 Electrolytes Based on Ionic Liquids Ionic liquids (ILs) have been examined in the preparation of electrolytes for Mg batteries because, in addition to being endowed with high thermal and electrochemical stability, they exhibit negligible vapor pressure and are non-flammable.150-151 However, it was determined that magnesium electrodes are generally reactive toward imidazolium-based ILs.152 Magnesium appears to develop a blocking passivation layer impervious to the transport of ions in the few ILs that are electrochemically stable. These observations dampened earlier interest in EMImCl/AlCl3 melts (EMImCl: 1-ethyl-3-methylimidazolium chloride) partially neutralized with conventional MgCl2 153-154 which were carefully studied due to their intriguing Lewis acid-base properties defined in terms of the ratio 𝑅 = 𝑛𝐴𝑙𝐶𝑙3 𝑛𝐸𝑀𝐼𝑚𝐶𝑙⁄ . A melt is classified as basic, neutral, or acidic, if R is less, equal to, or more than 1, respectively. High-performance electrolytes based on EMImCl doped with AlCl3 and highly amorphous MgCl2 have recently been developed.138 These systems have the general formula: (EMImCl/(AlCl3)1.5)/(- 19 MgCl2)x, with molar ratio 𝑥 = 𝑛𝛿−𝑀𝑔𝐶𝑙2 𝑛𝐼𝐿⁄ , 0  x  0.2, and exhibit an initial capacity of 80 mAh/g with a steady-state voltage of 2.3 V when discharged in a coin cell set up. The structural characteristics and the properties of the resulting electrolyte are tied to the ion-ion interactions within the IL. This development seems very promising, however the key to improvement and success of these new materials is identification of the structure. Experimental IR and Raman spectra of the electrolyte systems with the general formula (EMImCl/(AlCl3)1.5)/(-MgCl2)x reveals peaks which are indicative of the formation of new species and structural changes over the pristine IL, i.e., EMImCl/(AlCl3)1.5. Changing the halide from chloride to the iodide also affects the electrolyte properties. EMIm-based ILs have been studied extensively.155-159 However, the effect of the Mg-halide salt on their structure and properties has not been explored. Quantum chemical calculations were undertaken previously to determine the structure of small clusters of molecules and assignment of the experimental vibrational spectra in various ILs including: EMImCl;160-162 EMImCl/AlCl3; 160 EMImBr;163 EMImBF4; 163-164 and EMImPF6. 163 Hence, molecular-level investigations were undertaken in this work to understand the ionic structure and experimentally observed properties of these EMIm- based electrolytes. 1.4 Research Overview The research undertaken and reported in this dissertation aims to understand various structural and dynamical aspects of a variety of electrolytes used in energy storage and conversion technologies. Several modeling techniques at different length scales were utilized to elucidate how the molecular structure determines the properties of these materials. The hydration of the four vanadium cations (V2+, V3+, VO2+, and VO2 +) present in VRFB electrolytes was determined via electronic structure calculations. The Gibbs energies of hydration 20 of the four vanadium cations were calculated using quasi-chemical theory. Electronic structure calculations determined the local interactions of the hydrated vanadium cations with sulfuric acid (H2SO4) and trifluoromethanesulfonic acid (CF3SO3H) employing a continuum solvation model to account for the solvent/liquid phase effects. A parameterization method based on ab initio calculations was developed to compute interaction parameters utilized in dissipative particle dynamics (DPD) simulations of hydrated Nafion membrane with absorbed vanadium cations. The hydrated morphology of Nafion was investigated to determine the effects of absorbed vanadium cations. Microstructure/morphology relationships in the AEM systems based on SEBS was explored at the meso-scale via DPD simulations employing the developed parameterization scheme. The SEBS triblock copolymer was first benchmarked and then hydrated AEM systems were simulated to probe the effect of the degree of hydration and anion type (OH− vs Cl−) on the morphology. A coarse-graining scheme was proposed for the SDAPP membrane electrolyte systems. This was followed by electronic structure calculations in order to provide the atomic scale information required for the parameterization of the DPD interactions (utilizing the method developed in this work). Two recent EMIm-based electrolyte materials for magnesium batteries were investigated using density functional theory (DFT). The EMImCl/(AlCl3)1.5 and EMImAlI4 ILs were studied at the proximity of MgX2 (X = Cl and I, respectively) to investigate the effect of the magnesium salt. The calculated vibrational frequencies were compared to the experimental IR and Raman spectra to assign the vibrational peaks and reveal the speciation of the electrolyte systems. The hydration and dissociation of ten acid molecules that are commonly utilized in electrolytes: CF3SO3H; FSO3H; CH3SO3H; H2SO4; HNO3; HF; H3PO2; H3PO3; H3PO4; and CF3SO2NHSO2CF3 21 were studied via the electronic structure calculations. The hydrated structures were further explored for proton transfer by performing potential energy surface scans in both gas phase and with a dielectric continuum model. These calculations aimed to provide a fundamental understanding of the relative acidities and to serve the purpose of designing novel electrolytes featuring one or more of these acids in either the solution or tethered to an ionomer. The theoretical background for the various techniques employed in this work are summarized in Chapter 2. A discussion of the results is presented in Chapters 3 – 7: Chapter 3, multiscale modeling of VRFB electrolytes with quantum chemical calculations, quasi-chemical theory, and DPD simulations; Chapter 4, morphological investigations of SEBS-based AEMs coarse-grained DPD simulations; Chapter 5, a meso-scale study of SDAPP membranes; Chapter 6, structural and vibrational frequency analysis of EMIm-based electrolytes; and Chapter 7, DFT calculations of the hydration and proton transfer of several acids. The summary of the significant findings and concluding remarks are given in Chapter 8. 22 Chapter 2 Computational Methods 2.1 Multiscale Modeling Any physical phenomena may be viewed from different scales and levels of complexity depending on the application and properties of interest. For example, in engineering problems where the macroscopic behavior of materials is important, continuum models are often used with little or no appeal to molecular scale information. In contrast, in quantum physics the materials behavior is often described by the atomic and electronic level processes.165 Figure 2.1 shows a classification of various computational methods in chemistry and materials science as functions of length and time scales: quantum mechanics; atomistic; mesoscopic; and macroscopic. Decreasing the length scale usually increases the computational expense and thereby decreases the time scale. Figure 2.1. Classification of several computational methods in chemistry and materials science as a function of length and time scales. The methods at larger length scales probe longer times but lack any small scale information. Length Scale T im e S ca le s ms μs ns ps pm nm μm mm m 23 Macroscopic scale methods (e.g. finite element calculations) assume that the materials of interest are continuous and ignore the microscopic and atomic resolution. The effects of the smaller scales may be implicitly introduced in these methods through constitutive equations that are often empirical and based on simple considerations such as symmetry, linearization, or Taylor series expansions.165 Empirical constitutive models are typically highly successful in engineering applications, however they are inadequate in cases where the system is complex and lack the ability to connect information between the microstructure and macroscale behavior. Dissipative particle dynamics (DPD), Brownian dynamics (BD), and coarse-grained molecular dynamics (CG-MD) are meso-scale simulation techniques that have been extensively employed48- 50,54,133,166 and improved167-172 in recent years due to their computational efficiency. In these methods, several atoms or molecules are replaced with a representing mesoscopic particle (or bead) and a mesoscopic force field assumed to define the interactions between the particles. The mesoscopic force fields are critical to the coarse-grained models. Figure 2.2 shows a schematic of a coarse-grained model for the PFSA ionomer, Nafion. Figure 2.2. Coarse-graining of Nafion, a PFSA ionomer. Different mesoscopic particles are defined to represent chemically distinct subunits of the polymer (taken from reference 173). 24 Classical molecular dynamics and Monte Carlo (MC) are atomic scale simulation techniques that consider full atomistic resolution. The atoms and molecules interact with each other based on force fields that are obtained either empirically or by fitting to quantum mechanical data.174- 175 However, these methods do not consider electrons and subatomic particles and therefore do not effectively describe electronic effects including chemical reactions, breaking and forming of covalent bonds, and to some extent even hydrogen bonding. The methods at the finest scale shown in Figure 2.1 rest on the principles of quantum mechanics to determine the electronic structure of the matter. First principles (i.e., ab initio) calculations are the most accurate and computationally expensive methods and therefore limited to very small systems consisting of hundreds of atoms or several molecules. The application of quantum mechanics to the problems in chemistry is referred to as quantum chemistry.176 In what follows, a brief description of quantum chemical calculations and some approximation techniques are presented. Next, a statistical mechanical approach that describes the thermodynamics of ion solvation phenomena, quasi-chemical theory, is presented. The statistical mechanical derivations and practical formulations required for the calculation of Gibbs energies are presented. The meso-scale simulation method, DPD, is then described and several parameterization methods for obtaining the coarse-grained interactions are reviewed. 2.2 Quantum Chemical Calculations The insufficiency of Newtonian mechanics in describing the behavior of very small particles such as electrons and atoms nuclei led to the development of quantum mechanics. This branch of physics considers the discreteness of energy, the wave-particle duality of light and matter, and spin of a particle that Newtonian mechanics fails to describe. In quantum mechanics, the state of a physical system at a given time is postulated to be described in the form of a function 25 of the particles’ coordinates and spin, namely the wave function Ψ.176 The state of a system changes with time and therefore the wave function is also time dependent. The time evolution of the system is described by the Schrödinger equation that predicts the future state of the system. 𝑖ℏ 𝜕 𝜕𝑡 Ψ = 𝐻̂Ψ (2.1) Here t is time, 𝑖 = √−1 is the imaginary unit, ℏ is the plank’s constant divided by 2𝜋, and 𝐻̂ is the Hamiltonian operator that describes the total energy of the system. If the system experiences no time-dependent external forces and the potential energy is not a function of time, Ψ may be written as Ψ = f(t)𝜓, i.e., the product of a function of time, f(t), and a function of particles’ coordinates and spin, 𝜓 . In these cases the Schrödinger equation reduces to a time-independent form for 𝜓. Such states are called stationary states and do not change with time. Static ab initio electronic structures are based on the solution of the time-independent Schrödinger equation for a physical system of interest. However, the exact solution to this equation is only available for a single-particle system. Hence, different approximation methods are developed and used. The time-independent Schrödinger equation is introduced along with approximate solution methods including Hartree-Fock (HF) and density functional theory (DFT).176-177 Later the basis sets, a brief introduction to continuum solvation models, and vibrational frequencies are explained. 2.2.1 Time-Independent Schrödinger Equation In quantum mechanics the stationary state of a physical system is described by wave function, a function of the particles’ coordinates and spin. Every measurable physical quantity is described by a Hermitian operator (i.e., 𝐴̂ = 𝐴̂†). The operator corresponding to the total energy of the system 26 is the Hamiltonian operator, 𝐻̂. It is also postulated that the measurement of a physical quantity gives an eigenvalue of its operator, and therefore the physical quantity can be calculated by solving an eigenvalue problem. The time-independent Schrödinger equation is an eigenvalue problem for the total energy of the system: 𝐻̂𝜓 = 𝐸𝜓 , (2.2) where 𝜓 is the space and spin wave function/eigenfunction, and 𝐸 is the energy/eigenvalue. The wave function itself does not have any physical interpretation. However the product of 𝜓 with its complex conjugate 𝜓∗, which is |𝜓∗𝜓| = |𝜓|2 is the probability of finding the system in the state 𝜓, provided that the normalization condition: ∫|𝜓|2 𝑑𝜏 = 1 , (2.3) is satisfied. The 𝜏 is coordinates and spin of all particles and the ∫ 𝑑𝜏 integrates over all the coordinates and spin in which the system and its Hamiltonian are defined. The Hamiltonian is simply the sum of the kinetic and potential energies of all the particles in the system. 𝐻̂ = 𝑇̂𝑒 + 𝑇̂𝑛 + 𝑉̂𝑒𝑒 + 𝑉̂𝑛𝑛 + 𝑉̂𝑒𝑛 (2.4) Here 𝑇̂ and 𝑉̂ represent the operators of kinetic and potential energies, respectively and the subscripts e and n refer to electron and nucleus, respectively. The operator 𝑉̂𝑒𝑒 represents the potential energy due to the electron-electron interactions and similarly 𝑉̂𝑛𝑛 and 𝑉̂𝑛𝑒 are the operators corresponding to the nucleus-nucleus and nucleus-electron interactions. The expanded form of the Hamiltonian operator in atomic units for 𝑁𝑒 electrons and 𝑁𝑛 nuclei is: 27 𝐻̂ = − 1 2 ∑∇𝑖 2 𝑁𝑒 𝑖 −∑ 1 2𝑀𝑘 ∇𝑘 2 𝑁𝑛 𝑘 +∑∑ 1 𝑟𝑖𝑗 𝑁𝑒 𝑗>𝑖 𝑁𝑒 𝑖 +∑∑ 𝑍𝑘𝑍𝑙 𝑟𝑘𝑙 𝑁𝑛 𝑙>𝑘 𝑁𝑛 𝑘 −∑∑ 𝑍𝑘 𝑟𝑖𝑘 𝑁𝑛 𝑘 𝑁𝑒 𝑖 (2.5) where i and j denote electrons, k and l denote nuclei, 𝑀𝑘 is the ratio of the mass of the kth nucleus to the mass of an electron, ∇2 is the Laplacian operator, Z is the atomic number of the nucleus, and 𝑟𝑖𝑗 = |𝑟𝑖 − 𝑟𝑗| and 𝑟𝑘𝑙 = |𝑟𝑘 − 𝑟𝑙|. The mass of the nuclei is far greater than that of electrons, implying that the time scale of nuclear motion is significantly smaller than the electrons and hence their motion may be decoupled. Thus, in the study of the motion of electrons, the nuclei can be considered at fixed positions. This is the Born-Oppenheimer (BO) approximation and permits the neglecting of the kinetic energy of the nuclei in electronic structure calculations. Because the position of the nuclei are no longer variables, the nuclear-nuclear repulsion becomes a constant which solely depends on the positions of the nuclei. Hence, the Hamiltonian operator becomes: 𝐻̂ = 𝐻̂𝑒𝑙𝑒𝑐 + 𝑉𝑛𝑛 (2.6) where, 𝐻̂𝑒𝑙𝑒𝑐 is the electronic Hamiltonian. This may be further simplified by introducing the compact notation for the one- and two- electron terms: 𝐻̂𝑒𝑙𝑒𝑐 =∑ℎ𝑖 𝑁𝑒 𝑖 +∑∑𝑔𝑖𝑗 𝑁𝑒 𝑗>𝑖 𝑁𝑒 𝑖 . (2.7) Here, ℎ𝑖 = − 1 2 ∇𝑖 2 − ∑ 𝑍𝑘 𝑟𝑖𝑘 𝑁𝑛 𝑘 is the one-electron operator corresponding to the kinetic energy and coulombic attraction of electron i in the field of all the nuclei, and 𝑔𝑖𝑗 = 1 𝑟𝑖𝑗 is the two- electron operator describing the electron-electron repulsion. The electronic version of the Schrödinger equation may be written as follows: 𝐻̂𝑒𝑙𝑒𝑐𝜓𝑒𝑙𝑒𝑐 = 𝐸𝑒𝑙𝑒𝑐𝜓𝑒𝑙𝑒𝑐 , (2.8) 28 where 𝜓𝑒𝑙𝑒𝑐 and 𝐸𝑒𝑙𝑒𝑐 are the electronic wave function and electronic energy, respectively. 2.2.2 Hartree-Fock Method The difficulty in evaluating an analytic solution of the time-independent Schrödinger equation is the coulombic electron-electron repulsion term which prevents independent treatment of the electrons. To overcome the problem associated with electron correlation it is assumed that each electron moves in an orbital within the average field generated by all the other electrons. This approximation178 leads to a simple form of the many-electron wave function. The N-electron wave function is expressed as the product of N independent single particle wave functions. However, the electronic wave function must be antisymmetric in order to satisfy the Pauli Exclusion Principle, which states that each electron in a system has a unique quantum state, i.e., the probability of finding two electrons in the same state is zero. To satisfy the antisymmetric requirement the wave function is expressed as a Slater determinant (SD): 𝜓𝑆𝐷(𝑥1, 𝑥2, … , 𝑥𝑁) = 1 √𝑁! | 𝜒1(𝑥1) 𝜒2(𝑥1) ⋯ 𝜒𝑁(𝑥1) 𝜒1(𝑥2) 𝜒2(𝑥2) ⋯ 𝜒𝑁(𝑥2) ⋮ 𝜒1(𝑥𝑁) ⋮ 𝜒2(𝑥𝑁) ⋱ ⋯ ⋮ 𝜒𝑁(𝑥𝑁) | , (2.9) where 𝑥𝑖 is the space-spin coordinate for the ith electron, 𝑥𝑖 ≡ [𝑟𝑖 , 𝜔𝑖] (𝑟𝑖 and 𝜔𝑖 are the spatial and spin coordinates, respectively) and 𝜒𝑖 is the spin orbital. The spin orbitals must be orthonormal, i.e., ⟨𝜒𝑖|𝜒𝑗⟩ = 𝛿𝑖𝑗 and therefore the prefactor 1 √𝑁!⁄ ensures the normalization of the wave function. Using the SD and Equation (2.7) the electronic HF energy of the system may be calculated as: 𝐸𝑒𝑙𝑒𝑐 𝐻𝐹 = ⟨𝜓𝑆𝐷|𝐻̂𝑒𝑙𝑒𝑐|𝜓 𝑆𝐷⟩ = ⟨𝜓𝑆𝐷| ∑ ℎ̂𝑖𝑖 |𝜓𝑆𝐷⟩ + ⟨𝜓𝑆𝐷| ∑ ∑ 𝑔̂𝑖𝑗𝑗>𝑖𝑖 |𝜓𝑆𝐷⟩ . (2.10) By substituting a simpler form of the SD, using the orthonormality of the spin orbitals, and knowing the electronic HF energy can be written as the sum of one- and two- electron terms as: 29 𝐸𝑒𝑙𝑒𝑐 𝐻𝐹 =∑ℎ𝑖 𝑁 𝑖 + 1 2 ∑∑(𝐽𝑖𝑗 − 𝐾𝑖𝑗) 𝑁 𝑗 𝑁 𝑖 , (2.11) where the one-electron term ℎ𝑖 is the expectation value of its operator, i.e., ⟨𝜓𝑆𝐷|ℎ̂𝑖|𝜓 𝑆𝐷⟩ = ℎ𝑖, and corresponds to the kinetic energy and coulombic attraction of electron i in the field of all the nuclei. The term 𝐽𝑖𝑗 is called the Coulomb integral as it represents the coulombic repulsion between the two i and j electrons (𝐽𝑖𝑗 = ⟨𝜒𝑖(𝑥𝑖)𝜒𝑗(𝑥𝑗)|𝑔̂𝑖𝑗|𝜒𝑖(𝑥𝑖)𝜒𝑗(𝑥𝑗)⟩) and 𝐾𝑖𝑗 is the exchange integral and unlike the Coulomb integral has no physical analog (𝐾𝑖𝑗 = ⟨𝜒𝑖(𝑥𝑖)𝜒𝑗(𝑥𝑗)|𝑔̂𝑖𝑗|𝜒𝑗(𝑥𝑖)𝜒𝑖(𝑥𝑗)⟩). The exchange integral is the result of the antisymmetric requirement of the wave function. If the spin orbitals and consequently the SD (the wave function) are known the energy of the system can be calculated accordingly. Unfortunately, the mathematical form of the wave function is unknown in almost all cases. The practice is to apply a variational method and to make use of a trial wave function to calculate the energy of the system. This trial wave function is then optimized until the energy of the system reaches a minimum. The energy obtained via the trial wave function is always greater than the ground state energy unless the wave function is exact. The HF method applies the variational principle to determine the spin orbitals that minimize the energy (the expectation value of the Hamiltonian). The method of Lagrange multipliers under the constraint of orthonormal orbitals is applied to the energy expression in Equation (2.11). This leads to the HF equation for an arbitrary electron: 𝐹̂(1)𝜒𝑖(1) = 𝜀𝑖𝜒𝑖(1) , (2.12) where 𝐹̂ is the Fock operator defined as: 𝐹̂(1) = ℎ̂1 +∑(𝐽𝑗(1) − 𝐾𝑗(1)) 𝑁 𝑗 , (2.13) 30 which is the single electron operator describing the energy of an electron involved in the interaction with all other electrons and the nuclei. The first term of the Fock operator is associated with the kinetic energy of an electron and its attractive interaction with all the nuclei. The second term describes the electron-electron repulsion due to the presence of all other electrons. Both Coulomb and exchange operators depend on all the spin orbitals and hence the Fock operator itself. Thus, the Fock operator can be determined only if all other occupied orbitals are known. The consequence of this dependency is that the HF equation becomes a non-linear eigenvalue problem, which involves iterative procedures with an initial guess of the orbitals to solve. The iteration is usually continued until the change in the energy is less than a predetermined value. This iterative approach is referred to as the self-consistent field (SCF) method. The term 𝜀𝑖 in the Equation (2.12) can be interpreted as the energy of the ith molecular orbital, 𝜒𝑖(1). The molecular orbitals are usually expressed as a linear combination of atomic orbitals (LCAO) that are expanded as: 𝜒𝑖 =∑𝐶𝑎𝑖𝜑𝑎 𝑎 , (2.14) where 𝜑𝑎 is the atomic orbital of atom a (known) and 𝐶𝑎𝑖 is the expansion coefficient (unknown). The HF equations cast in matrix form under the LCAO approximation are the Hartree-Fock- Roothaan equations: 𝐹𝐶 = 𝑆𝐶𝜀 , (2.15) where 𝐹 and 𝑆 are the Fock and basis function overlap matrices, respectively ( 𝐹̂𝑎𝑏 = ⟨𝜑𝑎|𝐹̂(1)|𝜑𝑏⟩ and 𝑆𝑎𝑏 = ⟨𝜑𝑎|𝜑𝑏⟩). 𝐶 is the coefficient matrix and 𝜀 is a diagonal matrix of the orbital energies. A first set of coefficients along with the atomic basis functions are used as an initial guess to find the orbital energies in Equation (2.12) via the SCF method. 31 HF theory assumes a single-determinantal form for the wave function and neglects part of the correlation between the electrons. The electrons are subject to an average non-local potential arising from the other electrons, which can lead to a poor description of the electronic structure. Although the HF method provides a base for other more sophisticated post HF methods such as Møller-Plesset (MP) perturbation theory and configuration interaction (CI), but its accuracy is questionable. Post HF methods generally provide a better description of the electron correlation and consequently a more accurate description of the system, however they are considerably more computationally expensive. 2.2.3 Density Functional Theory DFT has been widely used in recent years to overcome the electron correlation deficiency and computational cost of the HF and post HF methods. In this approach the molecular electronic energy is obtained from the molecular electron probability density rather than the molecular wave function. This significantly decreases the number of variables, since the molecular electron probability density is a function of only 3 spatial coordinates for all the electrons, whereas the molecular wave function depends on 3N spatial and N spin coordinates. The central idea to DFT was first introduced by Hohenberg and Kohn in 1964.179 They developed two theorems in which the electronic energy could be obtained exactly from the electron density. In the first theorem, called the existence theorem, they proved that for a system of electrons in an external potential with a non-degenerate ground state, the ground state electron probability density has a one to one relationship with the ground state molecular energy, wave function, and all other molecular electronic properties. Thus, the ground state electronic energy, 𝐸0, can be written as a function of the ground state electron probability density, 𝜌0, as: 𝐸0 = 𝐸𝜐[𝜌0] = 𝑇𝑒[𝜌0] + 𝑉𝑒𝑒[𝜌0] + 𝑉𝑛𝑒[𝜌0] (2.16) 32 where the square brackets denote a functional relation and the subscript 𝜐 emphasizes the dependence of 𝐸0 on the external potential. The 𝑇𝑒 and 𝑉𝑒𝑒 are the ground state kinetic energy and electron-electron repulsion energy functionals, respectively. The 𝑉𝑛𝑒 functional is associated with the attraction between the electrons and the nuclei and is only a function of the coordinates of the electrons (the nuclei are fixed due to the BO approximation). This potential is considered external in DFT, since it is produced by the nuclei external to the system of electrons. If the electron density is known the 𝑉𝑛𝑒[𝜌0] can be determined, however the other two functionals, 𝑇𝑒[𝜌0] and 𝑉𝑒𝑒[𝜌0], are unknown. The second Hohenberg and Kohn theorem, the variational theorem, states that the exact ground state electron density should minimize the energy functional which is in fact the application of the variational principle (similar to the HF method). Thus, 𝐸𝜐[𝜌𝑡𝑟𝑖𝑎𝑙] ≥ 𝐸𝜐[𝜌0] (2.17) and the trial density must satisfy the ∫𝜌𝑡𝑟𝑖𝑎𝑙(𝑟)𝑑(𝑟) = 𝑁 and 𝜌𝑡𝑟𝑖𝑎𝑙(𝑟) ≥ 0 conditions. Although the Hohenberg-Kohn theorems provide the possibility to calculate all the ground state molecular properties from 𝜌0, they neither explain how to determine 𝐸0 without knowing the kinetic and repulsion potential energy functionals nor how to find 𝜌0 without first knowing the form of the wave function. Kohn and Sham devised a practical method to utilize the Hohenberg-Kohn formalism.180 They considered a hypothetical system, denoted by subscript s, consisted of N non-interacting electrons that each experience the same external potential energy function, where this potential is such as to make the ground state electron density of the hypothetical system equal to the ground state electron density of the real system, i.e. 𝜌𝑠 = 𝜌0 (the density subscript will be dropped for simplicity in what 33 follows). To facilitate analysis of the unknown functionals in Equation (2.16), the energy functional may be rewritten as: 𝐸𝜐[𝜌] = ∫𝜌(𝑟)𝜐(𝑟)𝑑𝑟 + 𝑇𝑠[𝜌] + 1 2 ∬ 𝜌(𝑟1)𝜌(𝑟1) 𝑟12 𝑑𝑟1𝑑𝑟2 + 𝐸𝑥𝑐[𝜌] , (2.18) where the first term is the external attraction potential energy, 𝑇𝑠 is the kinetic energy functional of the non-interacting system, the double integral term is the self-repulsion of a classical charge distribution, and the last is the exchange-correlation energy functional, 𝐸𝑥𝑐[𝜌]. The exchange- correlation functional includes not only the effects of exchange and correlation embedded in 𝑉𝑒𝑒[𝜌], but also the correction for the classical self-interaction energy (the double integral term) and for the difference in kinetic energy between the fictitious non-interacting system and the real one (Δ𝑇[𝜌]). Determining the electron density for the hypothetical non-interacting system is easy because there is no electron-electron repulsion term and the electronic Hamiltonian can be easily expanded as the sum of one-electron operators and the eigenfunction is the SD of the individual one electron eigenfunctions. Rewriting the energy expression as: 𝐸[𝜌] = −∑⟨𝜒𝑖| ∑ 𝑍𝑘 𝑟𝑖𝑘 𝑁𝑛 𝑘=1 |𝜒𝑖⟩ 𝑁 𝑖=1 −∑⟨𝜒𝑖| 1 2 ∇𝑖 2|𝜒𝑖⟩ 𝑁 𝑖=1 +∑⟨𝜒𝑖| 1 2∫ 𝜌(𝑟𝑗) 𝑟𝑖𝑗 𝑑𝑟𝑗 |𝜒𝑖⟩ 𝑁 𝑖=1 + 𝐸𝑥𝑐[𝜌] , (2.19) where the density is 𝜌 = ∑ ⟨𝜒𝑖|𝜒𝑖⟩ 𝑁 𝑖=1 . In order to find 𝜒𝑖 a similar approach to the HF method is used. The method of Lagrange multipliers is applied to minimize the energy expression in this equation and results in the central KS DFT equations: ℎ̂𝑖 𝐾𝑆𝜒𝑖 = 𝜀𝑖𝜒𝑖 , (2.20) where ℎ̂𝑖 𝐾𝑆 is the one-electron KS Hamiltonian: 34 ℎ̂𝑖 𝐾𝑆 = − 1 2 ∇𝑖 2 −∑ 𝑍𝑘 𝑟𝑖𝑘 𝑁𝑛 𝑘 +∫ 𝜌(𝑟𝑗) 𝑟𝑖𝑗 𝑑𝑟𝑗 + 𝑉𝑥𝑐 . (2.21) 𝑉𝑥𝑐 is the functional derivative: 𝑉𝑥𝑐 = 𝛿𝐸𝑥𝑐 𝛿𝜌 . (2.22) The KS equations are solved iteratively, similar to t