ATOMIC SIMULATION ON CHEMICAL-MECHANICAL COUPLED DEFORMATIONS IN COMPLEX NANO STRUCTURES By Jialin Liu A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Material Science and Engineering – Doctor of Philosophy 2019 ATOMIC SIMULATION ON CHEMICAL-MECHANICAL COUPLED DEFORMATIONS IN COMPLEX NANO STRUCTURES ABSTRACT By Jialin Liu Nano-structured materials often exhibit very different mechanical properties comparing with their bulk counterpart and are more sensitive and active to chemical interactions with the environments due to the large surface-to-volume ratio. In this thesis, predictive modeling techniques including density functional theory (DFT) and reactive molecular dynamics method (MD) are designed and applied to understand the deformation mechanisms of complex nano-structured material and describe chemical-mechanical coupled interactions. Three technologically important materials are investigated, to understanding the high strain rate toughening mechanism in nacre, predicting the formation and fracture of aluminum oxide bifilms in aluminum castings, and revealing the lithium growth morphology as a function of oxygen partial pressure. For nacre, its hierarchical structure and toughening mechanisms have inspired many materials developments. Recently, a new toughening mechanism, deformation twins was observed in nacre after dynamic loading (103 s-1). The deformation twinning tendency and the competition between fracture and deformation twinning were revealed by DFT calculations. We discovered that the ratio of the unstable and the stable stacking fault energy in aragonite is hitherto the highest in a broad range of metallic and oxide materials and the bonding nature for this high ratio is explained. Both aluminum and lithium have high oxygen affinity. Their interaction with the oxygen environment affects the mechanical properties and vice versa. During casting of aluminum, it has long been proposed that the entrapped alumina “bifilms” are detrimental to the fatigue properties of the cast product. However, its properties have never been measured due to experimental limitations. Therefore, a ReaxFF based MD protocol was designed to simulate aging, folding, and fracture of oxide bifilms. The predicted fracture energy, fracture location, and differences between old and young oxides are explained a series of experimental observations. To illustrate the Li- growth mechanism in a solid-state-battery testing platform, we modeled the morphology of Li nano-structure growth in oxygen environment via ReaxFF-based MD. The simulation revealed that the competition of the Li growth rate and oxidation rate leads to the sphere-nanowire-sphere morphology transition with increasing oxygen partial pressure. Understanding the impact of chemical reaction on Li dendrite growth mechanisms and morphology evolution provided insights on the formation of the solid electrolyte interface (SEI) layer in a Li-ion battery. Finally, a shortcoming of the current charge transfer scheme (qEq) used in the ReaxFF MD simulation is discussed. It is demonstrated that qEq method will lead to overductile ionic materials in the MD simulation. A new Force field method and new parameters are proposed to mitigate this problem. ACKNOWLEDGEMENTS Since the year 2014, I joined Prof. Yue Qi’s group with little knowledge of research and my ignorance. As a careless person, I constantly made many mistakes in the business that requires great amount of carefulness. Although the motto “details make a difference” lies in my brain for several years, I learned it again in the hard way by realizing none of my jobs can run successfully without correcting some parameters at the first time… With this valuable lesson learned on the cost of many wasted CPU time, I want to express my greatest gratitude to my advisor Prof. Yue Qi for all her great amount of patience and kind guidance in these five years. In the five years, it is my horror to have many opportunities for collaborating with other excellent research groups which helped me to develop better communication skills and provided me more insight for understanding science from other perspectives. I’d like to thank Prof. Matthew Hirn for helping me on the machine learning project, Prof. Qigui Wang from GM for the aluminum casting project, Prof. Martin Crimp for the Ni 9R project, Prof. Metin Aktulga for forcefield fitting project, and Prof. Wei Lai for part of the SPG project. Also, my study experience in MSU is made better because of the generous help from many faculties, including my committee members, Prof. Andre Lee, Prof. Thomas Bieler and Prof. Lixin Dong. I have been enjoying taking or auditing many classes in MSU which not only helped with my research work but also extended my vision to other related fields. Specifically I want to sincerely thank Prof. Kalyanmoy Deb for optimization method, Prof. Thomas Bieler for dislocation, Prof. Yue Qi for computational modeling, Prof. Jiayu Zhou for machine learning, Prof. Martin Crimp for TEM, Prof. Matthew Hirn for math foundation in data science, and all other classes that I have attended. iv Lab is the place where I spend most of the time in. My lab mates provided me great support for many intelligent discussions on both struggled research and life. I would like to thank Dr Sungyup Kim, Dr. Christine James, Dr. Tridip Das, Dr. Kwangjin Kim for helping me getting familiar with our lab and the MSU campus when I first arrived in MSU. Thank Dr. Thanaphong Phongprecha (Joe) for being both my lab mate and roommate, and great amount of trust. Thank my neighbor Aritra Chakraborty for helping solving all kinds of strange research problems and life problems. Thank Yuxiao Lin, and Hongkang Tian for walking with me through the five years of PhD life. Thank Dr. Yunsong Li, Dr. Michael Swift and Dr. Chi-ta Yang for giving me enormous help on my research work. The research work is, most of the time, enjoyable but sometimes boring. I feel lucky to have a bunch of friends share same hobbit with me. Thank the Chinese student community in Engineering Building for being my lunch pal and staying together to survive the first year culture shock including Yifeng Tian, Xing Lu, Zeyang Yu, Pei Chen, Yuxi Ma, Mingmin Wang. Dr. Quan Zhou and Dr. Yuanchao Liu. Thank all my skiing and snowboarding friends in Michigan snow team (MIST) for filling my world with shredding mountains and fresh powders. Thank all my climbing mates for spending time with me in Planet Rock and helping me defeating the fear of height, falling off, and myself. Thank my training partner Yaojie Liu for kicking me into the world of Brazilian Jiu Jitsu and my coach Matthew Linsemier in Majic BJJ for not only teaching me great techniques but also making the gym like a cozy home. Thank Aritra Chakraborty for dragging me to the world of CrossFit and my coaches Anthony Nye and Katie Trinklein in Greater Lansing CrossFit (GLC) for professional training advices and inspiring whips in the WOD. Thank my coaches Jimmy Wang and Huping Tian from CUSA Gyrfalcon Combat Academy for organizing many joyful wild boxing training sessions and wild chicken wing sessions. v Last but not least, I’d like to thank my parents for the continuous and unconditional support and love in the past 28.03 years. 余自甲午年中拜于齐先生门下，承蒙齐先生传道受业解惑至今已四载有余。遥想当 年，吾资质平平且懒惰愚笨，幸蒙先生不弃，收入门下。之后求学之路虽有艰难险 阻常常荆棘丛生，但幸有良师不厌其烦谆谆教诲，得以稳扎稳打屡屡逢凶化吉。为 求因材施教，先生特拿下新算法研发之项目以契合吾兴趣之所在，有师如此，实为 学生之幸。如今终得侥幸完成学业，日后必谨遵师嘱，刻苦钻研，不求功名显赫， 但愿砥砺前行，不忘初心。 vi TABLE OF CONTENTS CHAPTER 1. Background ....................................................................................................... 1 1.1. Complex deformation mechanisms in nano-structured materials ................................... 2 1.2. Multiscale modeling techniques on deformation mechanisms ..................................... 12 Individual Deformation mechanism study using DFT.......................................... 13 1.2.1. 1.2.2. Deformation process and mechanisms studied using the traditional forcefield based molecular dynamics ............................................................................................................... 16 1.2.3. Deformation mechanism studied using FEM ....................................................... 22 1.3. Motivation and Dissertation Outline ............................................................................. 23 2.1. 2.2. 2.2.1. 2.2.2. 2.2.3. CHAPTER 2. Chemical bonding origin of deformation twinning formation in the biomineral aragonite ................................................................................................................... 29 Summary ....................................................................................................................... 29 Introduction ................................................................................................................... 29 Hierarchy structure of nacre .................................................................................. 29 Toughening mechanisms of nacre under quasi-static loading .............................. 33 Experimental observations under dynamic loading .............................................. 36 2.3. Model Construction and Computational Methods ........................................................ 40 Compute the GSFE and Twinning Tendency ....................................................... 40 A short introduction to DFT method .................................................................... 42 Construction of nacre slab and computation details ............................................. 47 Results and Discussion ................................................................................................. 49 Surface Energy of Aragonite CaCO3 .................................................................... 49 Generalized SFE and Fracture in CaCO3 .............................................................. 51 Compare the twinning tendency of CaCO3 with other materials .......................... 53 The Atomic Structure and Bonding for the high USFE/SFE ratio ....................... 56 Conclusion .................................................................................................................... 59 2.4.1. 2.4.2. 2.4.3. 2.4.4. 2.3.1. 2.3.2. 2.3.3. 2.4. 2.5. 3.3.1. 3.3.2. 3.3.3. CHAPTER 3. Coupling the Formation Chemistry and the Mechanical Fracture of Oxide Bifilms in Cast Aluminum............................................................................................... 60 Summary ....................................................................................................................... 60 3.1. 3.2. Introduction ................................................................................................................... 61 3.3. Model Construction and Computational Methods ........................................................ 65 Forcefield based MD methods .............................................................................. 65 Simulation of the forming process of oxide bifilms ............................................. 70 Simulation of the fracture process of the aluminum oxide bifilms ....................... 74 3.4. Validation of ReaxFF .................................................................................................... 74 Scaled temperature ................................................................................................ 74 Density and Elastic Modulus of the oxides........................................................... 76 Results and Discussion ................................................................................................. 80 Interface analysis in oxide bifilms ........................................................................ 80 Deformation of the aluminum oxide bifilms under tension .................................. 85 Scale bridging to microscale for cohesive zone model ......................................... 89 3.5.1. 3.5.2. 3.5.3. 3.4.1. 3.4.2. 3.5. vii 3.6. Conclusion .................................................................................................................... 92 4.1. 4.2. 4.2.1. 4.2.2. CHAPTER 4. Oxidation modulated Lithium nano-structure growth morphology change during extrusion ............................................................................................................. 94 Summary ....................................................................................................................... 94 Introduction and Experimental Observation ................................................................. 95 The Li metal and Li2O in battery applications ...................................................... 95 Experimental observations .................................................................................... 97 Computational details ................................................................................................... 99 Oxidation of Li and Al slab .................................................................................. 99 Extrusion of Li nanowire .................................................................................... 101 Different conditions in MD simulation ............................................................... 102 Results and discussion ................................................................................................ 103 4.4.1. Validation of ReaxFF reactive forcefield ........................................................... 103 4.4.2. MD simulation of the oxidation process on Li and Al........................................ 104 Li growth morphology in oxygen environment .................................................. 109 4.4.3. Competition between oxidation rates and growth rates ...................................... 113 4.4.4. Conclusion .................................................................................................................. 115 4.3.1. 4.3.2. 4.3.3. CHAPTER 5. 4.3. 4.4. 4.5. 5.1. 5.2. 5.3. 5.4. 5.5. Improvement of charge transfer scheme in ReaxFF ................................. 116 Summary ..................................................................................................................... 116 Forcefield and its limitation on fracture simulation .................................................... 116 Computational approach ............................................................................................. 122 Results and discussion ................................................................................................ 123 Relaxation of the slab model............................................................................... 123 Deformation of the notched Li2O slab model .................................................... 125 Future work ................................................................................................................. 127 5.4.1. 5.4.2. CHAPTER 6. Summary and outlook .................................................................................. 128 viii LIST OF TABLES Table 2.1 Comparison of USFE, SFE UTFE and TFE values (mJ/m2) in metals, Si, perovskite oxides, calcite and aragonite, which has the highest (cid:1)(cid:2)(cid:3)(cid:4)(cid:5)/(cid:1)(cid:3)(cid:4)(cid:5) ratio .................................... 55 Table 3.1 Simulation details for monolayer and oxide bifilm structures shown in Figure 3.4 ..... 72 Table 4.1 The expected Li growth morphology at different charging rate and oxygen pressure. ..................................................................................................................................................... 114 ix LIST OF FIGURES Figure 1.1 A mechanism map for deformation behavior showing the nominal changes in the underlying mechanisms of plasticity at different grain sizes [48]. ................................................. 3 Figure 1.2 Schematic of the bivalve molluscan shell anatomy. It is composed of periostracum (a prismatic layer) and nacre. The liquid-filled interlamellar space exists between the mineralized shell and the mantle part of the soft body of the organism. The schematic also illustrates successive magnification of the brick and mortar structure of nacre [52]. ............................................................ 4 Figure 1.3 Photomicrographs from fracture surface of superplastically deformed AA5083 specimens tested in (a) air and (b) vacuum. In air the ligment is formed on fracture surface while in vacuum it showed a clean surface. .............................................................................................. 6 Figure 1.4 Computational and experimental techniques for a variety of length and time scales [68] ......................................................................................................................................................... 8 Figure 1.5 Comparison of the heating effect on (a,b) hydrogen-free and (c,d) hydrogenated pillars using ETEM. After annealing, the dislocations escaped from the surface for both pillars. But the hydrogenated pillars also formed a cavity between metal and surface oxides on the top of this pillar. ....................................................................................................................................................... 10 Figure 1.6 The multiscale simulation techniques for materials modeling at different length scale and time scale. Methods applied to the larger system usually need to compromise on calculation accuracy ........................................................................................................................................ 13 Figure 1.7 Typical terms contained forcefield for organic materials. The total energy consists of bond stretching, bond angle bending, bond twisting, van der Waals interactions, and Coulomb interactions. [105] ......................................................................................................................... 16 Figure 1.8 Proposed deformation mechanism map. σ/σ∞ is the normalized stress; r0/(cid:10) is the inverse of normalized grain size. [129] ......................................................................................... 19 Figure 1.9 Stacking fault formation in deformed aluminum nanowire in (a, b) vacuum and (c, d) O2 environment. Color bar indicates the central symmetry parameters. More dislocation nucleation events are observed in the O2 environment. .................................................................................. 21 Figure 2.1 Top: the hierarchical structure of nacre[52]. Bottom: The first order structure is a layered structure. The second order structure is platelet that’s further divided into smaller sub- plates [159].................................................................................................................................... 31 Figure 2.2 The TEM imagine of platelet structure in nacre shell [160] ....................................... 32 Figure 2.3 A schematic illustration of nacre shell structure including the platelet, biopolymer and mineral bridges [162] .................................................................................................................... 33 x Figure 2.4 Illustration of feature of platelet protruding in interlock mechanism [53] The platelets are locked at the circles to prevent relative displacement and rotation ........................................ 34 Figure 2.5 The short/long molecule model [54] (a) Schematic drawing of the biopolymer chain between two platelets. The long molecule chains on the top can be folded to modules and become short chain on the bottom. (b) The response of a short/long molecule chain under extension. The elongation can be increased while maintaining the flow stress by opening up modes on the molecule chain. ............................................................................................................................. 35 Figure 2.6 Experimental observation after deformation of nacre (a) Stress curve of nacre shell under different deformation rate and the HRTEM imagine of specimens after deformation. After quasi-static deformation, only (b) grain rotation is observed. But after dynamic deformation, (c) twinning and (d) partial dislocation emission can be observed. [146] ......................................... 37 Figure 2.7 Nacre’s extraordinary strengthening and toughening and its multi-functional deformation mechanisms. (a) The normalized toughness performance of different materials under high strain rate (~103) and low strain rate (~10-3). Nacre shows dramatic toughness increase over metals, polymers and even traditional composites under high deformation rate. (b) The deformation mechanisms that can be activated at the atomic level. When the deformation rate is high, only the mechanisms illustrated in grey platelets can be activated. .................................... 39 Figure 2.8 The atomic structure of deformation twinning (a) Upon impact along nacre’s cross section, deformation twinning/stacking fault takes over the deformation process, rather than the full-dislocation slipping[146]. (b) The atomic-scale measurement of white box area in (a) clearly demonstrates the twinning relationship between (110) and (110) planes. (c) The atomic 8-layer slab model is built according to (b) and is shown in two projections ........................................... 48 Figure 2.9 The surface energies of different layers and their relative errors comparing with 6-layer slab model. .................................................................................................................................... 50 Figure 2.10 The computed generalized stacking fault energy (GSFE) landscape and the fracture energy of aragonite. The USFE is higher than the fracture energy during rigid sliding, but drops below the fracture energy when layer expansion is allowed. Moreover, the USFE/SFE ratio is the highest in a broad range of metallic and oxide materials (as shown in Table 2.1). ...................... 52 Figure 2.11 Structural origin of different USFE/SFE ratio and the relaxation mechanism of stacking fault. Perfect slab, USFE and SFE structures of (a) aragonite, (c) Cu and (d) SrTiO3 show the lattice structure and Burger’s vector combination leads to high USFE/SFE ratio. The orange line is indicating the slipping plane. (b)Detailed analysis of one diamond shape crystal cell 2- in the short diagonal and maintain the same mass center indicates Ca2+ prefers to bond with CO3 in each CaCO3 column. In the fault plane, the equal tendency of jumping to two sites (crosses) 2- groups and dots in square and round lead to equalization of CO3-Ca bond. The triangles are CO3 are Ca2+. Different transparency indicates different depth ........................................................... 58 Figure 3.1 Schematic illustration of oxide bifilms formed during aluminum casting and their evolution with time and temperature change. ............................................................................... 65 xi Figure 3.2 Schematic of the ReaxFF iteration. The non-bonded interactions are on the left and the covalent/bonded interactions are on the right [131] ..................................................................... 68 Figure 3.3 The bond order of a carbon-carbon interaction with respect to distance. [131] .......... 69 Figure 3.4 The procedure to build aluminum oxide bifilm structures with different formation history ........................................................................................................................................... 71 Figure 3.5 Build structure for two-phase simulation. The initial temperature of 200K is applied to show the moving of interface toward crystal side ........................................................................ 76 Figure 3.6 The ReaxFF predicted room temperature Young's modulus of Al oxide with different oxygen content at different density ............................................................................................... 78 Figure 3.7 The predicted room temperature bulk modulus of bulk aluminum oxides at different phases in comparison with experimental values in the literature [233,259–270] ......................... 79 Figure 3.8 The charge on each atom of (a) the young oxides, (b) old oxides without extra heating, (c) old oxides and (d) –OH terminated old oxides with respect of their position along the surface normal direction. The corresponding atomic structures are shown .............................................. 81 Figure 3.9 (a) The Al-O bond density of (b) young oxides, (c) old oxides and (d) -OH terminated old oxides and their bifilm structures............................................................................................ 83 Figure 3.10 Deformation analysis of bifilms. (a) The fracture feature of young, old and –OH terminated old oxide bifilms. (b) The stress-strain curves of the three bifilms. ........................... 86 Figure 3.11 The fracture strengths of the three bifilms under different mesh size in finite element methods. ........................................................................................................................................ 90 Figure 4.1 (Left) A schematic of the experimental setup used to cycle the all-solid-state batteries in an operando SEM with the controlled pressure of O2. The grounded tungsten (W) tip is used to contact the top C-Al anode. A potential is applied to the bottom Pt-Ti current collector electrode to charge/discharge the electrochemical cell. (Right) Real-color optical image of C-Al anode in a representative device. .................................................................................................................... 96 Figure 4.2 Effect of oxygen pressure on Li plating morphology. SEM images of Li plated at (a) 5.7 × 10–7, (b) 5.7 × 10–6, and (c) 4.8 × 10–5 Pa residual oxygen pressure, respectively; 0.77 mA/cm2 current density. The SEM images are acquired at ≈ 0.04 mAh/cm2 capacity of deposition Li nanowire growth mechanism and critical thickness of the oxide layer .................................... 98 Figure 4.3 Schematic for slab model in MD simulation. The growth of Li nanowire is simulated by applying stress in oxygen environment. The oxidation rate is controlled by applying different temperature. ................................................................................................................................ 101 Figure 4.4 Surface morphology evolution of Li slab under oxidation at different time step. For each structure, the top view and the side view is showed. The red atoms indicate oxygen and purple atoms indicate Li. Only Li and O in the oxides are shown for a clear view. .............................. 105 xii Figure 4.5 Surface morphology evolution of Al slab under oxidation at different time step. For each structure, the top view and the side view is showed. The red atoms indicate oxygen and grey atoms indicate Al. Only Al and O in the oxides are shown for a clear view. ............................. 106 Figure 4.6 Number of oxygen atoms reacted with metal during the oxidation process. The Al slab shows more passivating effect than Li. ....................................................................................... 107 Figure 4.7 MD simulation snapshots of Li grown under (a) UHV, (b) oxygen atmosphere, and (c) oxygen atmosphere at an increased temperature to reach a higher oxidation rate after 50 ps for parts a and c and 55 ps for part b, respectively ........................................................................... 110 Figure 4.8 Schematics showing models of the Li nucleation and growth for the images presented in Figure 4.7. In the presence of O2, a lithium-oxide sheath is developed. Schematics are drawn not to scale. ................................................................................................................................. 111 Figure 5.1 Numerical results of charge on hydrogen atom for the dissociation of hydrogen fluoride. [153] ............................................................................................................................................ 121 Figure 5.2 The (a) as build structure of Li2O. Li atoms are in purple. Oxygen atoms are in red. ..................................................................................................................................................... 122 Figure 5.3 Relaxed Li2O slab structure using (a) qEq and (b) ACKS2 method. Li atoms are in purple. Oxygen atoms are in red. ................................................................................................ 123 Figure 5.4 Radial distribution function of two Li2O crystal structures and two relaxed slab structures ..................................................................................................................................... 124 Figure 5.5 Deformation of Li2O slab model using qEq method under (a) 1K and (b) 200K at different strains. Deformation with lower temperature showed more ductility during fracture. 125 Figure 5.6 The stress-strain curve of Li2O slab deformation process using the qEq method at 1K and 200K ..................................................................................................................................... 127 xiii CHAPTER 1. Background Materials with a feature length below 100 nm are classified as nano-structured materials [1]. Due to the large surface-to-volume ratio at the nanoscale and the high content of surfaces and interfaces, nano-structured materials exhibit unique mechanical properties that are very different from those with features sizes well above 100 nm. These properties include high fracture strength [2,3], modulus [4], and toughness [5,6]. For example, a Si nanowire with a diameter around 20 nm is a 1D nano-structured material widely used in electronic and thermoelectric devices [7,8] with high strength. The ultrathin 2D materials, such as graphene[9,10], MoS2[11,12], h-BN [13,14], MXene[15] etc, provide flexibility, strength, and possibly enticing electronic properties for innovative applications, such as spintronics and quantum information processing [16–18]. Materials found in nature often have hieratical structures based on nano-size building blocks. For example, the building block in nacre is an aragonite grain with a diameter of 20 – 45 nm surrounded by biopolymers such as chitin [19]. The ingenious structure of nacre is responsible for its high fracture toughness. A polycrystalline metal with a crystallite size less than 100 nm can be considered 3D nanostructured materials. They fill the gap between amorphous materials without any long-range order and conventional coarse-grained materials. The behavior of defects and mechanisms of plastic deformation in nanocrystalline metals and conventional polycrystalline metals are rather different. These unique mechanical properties of nano-structured materials have attracted much research interest because nano-structured materials can be used as building blocks of electronic, optical and nanoelectromechanical devices [20–27]. This new research calls for breakthroughs in novel deformation mechanisms that lead to high-performance materials. 1 1.1. Complex deformation mechanisms in nano-structured materials From the long history of metal processing, metallurgists have always been aware that reducing the grain size can significantly increase the mechanical performance of a metal such as hardness, toughness, and strength [28,29]. With the development of nanotechnology, many nanoscale materials exhibit dramatically different deformation behavior than their bulk counterparts. This phenomenon is often quantified by Hall-Petch relation [30–33]. For example, when the grain size of Cu drops from 20(cid:14)(cid:15) to 22 nm, its yield strength increases from 290 MPa to 480 MPa [34]. This can be captured by Equation 1.1: (cid:16)(cid:17)=(cid:16)(cid:19)+(cid:21)(cid:17)(cid:10)(cid:22)(cid:23)(cid:24) Equation 1.1 where (cid:16)(cid:17) is the yield stress, (cid:16)(cid:19) is the frictional stress, (cid:21)(cid:17) is a positive constant and (cid:10) is the grain size. The reason for the Hall-Petch relation is that a smaller grain size imposes a higher confinement effect on dislocation nucleation and movement, leading to harder and stronger nano- structured materials. However, when the grain size is smaller than a critical length (e.g. ~15 nm for Cu [34]), this mechanism breaks down and results into an inverse Hall-Petch relation. The nano-structured materials can become softer with decreasing grain size [31,35–37]. This is because when grains are too small, it is difficult to nucleate a dislocation. Due to the large grain-boundary- to-volume ratio in nanocrystalline metals, grain boundary diffusion and movement (including grain boundary sliding and grain boundary rotation) become the dominating deformation mechanism over dislocation interactions [38–41]. This deformation mechanism transition can also be triggered by the strain rate. For example, with nanoscale grain size, at appropriate high 2 temperature and slow strain rate, some aluminum and titanium alloys can show superplasticity with as much as 2000 percent deformation [42–45]. This is due to the plasticity enhanced grain boundary diffusivity [46,47]. The mechanism transition from Hall-Petch relation to inverse Hall Petch relation for nano-structured metals is illustrated in the schematic in Figure 1.1 [48]. In regions 1 and 2, the Hall-Petch relation holds because the full and partial dislocation generation and motion are the dominating deformation mechanisms. In region 3, it shows inverse Hall-Petch relation because the grain boundary dominates the plastic deformation. Figure 1.1 A mechanism map for deformation behavior showing the nominal changes in the underlying mechanisms of plasticity at different grain sizes [48]. 3D nano-structured inorganic materials can have even more complex hierarchical microstructures, enabling highly coordinated deformation mechanisms. Bio-materials such as nacre [49] and conch [50] have three to four levels of hierarchical structures as illustrated in Figure 1.2 (and Figure 2.1). 3 The first order structure of nacre is a layered structure with the thickness around 1 μm. Then it is divided into second-order structure which consisted of small platelets with a diameter of 1(cid:28)5 μm. The thickness of platelets is around 200 – 500 nm. Biopolymers such as chitin filled between platelets act as an adhesive. Furthermore, the platelets can be divided into a third-order structure, sub-platelets with pseudo-hexagonal platelets with a diameter of 30 – 80 nm. The third order structure was considered as the lowest level structure until Huang et al. [51] discovered the fourth order structure, nanoscale grains with a diameter of 20 – 45 nm. Also, their major component, aragonite CaCO3, has a lower crystal lattice symmetry and mixed ionic and covalent bonds, adding more complexity to the nanostructure. Figure 1.2 Schematic of the bivalve molluscan shell anatomy. It is composed of periostracum (a prismatic layer) and nacre. The liquid-filled interlamellar space exists between the mineralized shell and the mantle part of the soft body of the organism. The schematic also illustrates successive magnification of the brick and mortar structure of nacre [52]. 4 Although nacre is composed of 95 vol% aragonite CaCO3 and 5 vol% organic materials, its fracture toughness is much higher than that of CaCO3, because of its integrated hierarchical structure. Different deformation mechanisms corresponding to the hierarchical structures are discovered in nacre. At the platelets level, macroscopic deformation mechanisms, such as interlock mechanism [53] and the short/long molecule model [54], have been proposed to explain the extraordinary toughening properties of nacre. According to the interlock mechanism, the edges of each platelet are not smooth but have small kinks that can lock with kinks on another platelet to prevent the relative rotation. During deformation, these interlocks require extra energy to break. The short/long molecule model states that the biopolymer consists of short and long polymer chains. The short ones can provide high strength at smaller strain and absorb energy when they break. The long chains are folded into multiple modes that can be unfolded under stress to maintain the integrity of structure and strength at high strain. This mechanism has been used to explain the step-wise behavior during the stretching of bio-polymer. These combined features and deformation mechanisms lead to a remarkable combination of stiffness, high strength per weight, and high toughness per weight. The 1D and 2D nano-structured materials’ mechanical performance experience much greater influence by the environment, due to the large exposed surface-to-volume ratio. Since the most stable state of most elements is as an oxide, the surface of metals and alloys are often covered by an oxide layer when exposed to air. Early reviews reported that surface oxide tend to increase the yield stress of the metals and change the work hardening and fatigue strength [55]. The surface chemical reactions, such as oxidation, can also be coupled with the mechanical deformation process. Chang et al. reported the morphology difference between quasi-static stretching of fine grain Al-Mg AA5083 alloy in air and in vacuum [56]. As shown in Figure 1.3, the nanoscale 5 ligament can be seen on the fracture surface of the sample deformed in air while it showed a clean fracture surface in the vacuum. This work proved that the formation mechanism of the ligament is a mechanism of dynamic oxidation of alloy. Figure 1.3 Photomicrographs from fracture surface of superplastically deformed AA5083 specimens tested in (a) air and (b) vacuum. In air the ligment is formed on fracture surface while in vacuum it showed a clean surface. As the surface-to-volume ratio increases, the influence of the oxidizing environment on the mechanical properties becomes more dominant. This effect is especially significant for metals with high affinity to oxygen, such as Al, Mg, Li, and Ti. For example, even under an ultrahigh vacuum, e.g., at 10-8 Torr sec O2, an oxide layer can form on Al within seconds. In previous studies of mechanical properties of nano-wires, both in-situ transmission electron microscopy (TEM) and molecular dynamics (MD) simulations were often conducted in vacuum chambers, ignoring the environmental effect. But in practical application scenario, the chemical reaction between the material’s surface and the environment is often inevitable, especially for metals with high oxygen affinity such as Li and Al [57,58]. Due to the large surface-to-volume ratio in nano-structured 6 material, the chemical reactions occurring on the surface play more critical roles on the mechanical properties than the bulk part. It is now well accepted that the plasticity of metallic NWs deformed under vacuum follows a relationship of “smaller is stronger,” as its plastic deformation is controlled by the nucleation and escape of dislocations from the free surfaces. In contrast, Sen et al. reported oxidation induced softening and superplasticity in aluminum nanowires in the year 2014 via reactive MD simulations [59] and the results have been confirmed recently by environmental in situ TEM experiments in the year 2018 [60]. These new experiments and simulations have opened up many new questions on chemical-mechanical coupling in various nano-structures. One of the nano-structures that will be investigated in this thesis is the oxide-bi- films formed in cast aluminum [61]. The “environment” includes gas, such as air or oxygen, and liquid, such as water or even liquid electrolyte. The chemical reactions occurring on the surface involve oxidation and reduction reactions in more complicated operating environment. One particular environment for Li-metal is the electrolyte inside a Li-ion battery, containing LiPF6 salt dissolved in liquid ethylene carbonate and dimethyl carbonate. Li metal is not chemically stable in most organic electrolytes used today, just like metals are not chemically stable in air. The 10- to 100-nm thick passivation layer formed on the Li-metal electrode is called solid electrode interface (SEI). The major component of SEI is the decomposition product of electrolyte. The SEI ideally should allow the ion diffusion but block electron transportation, as electron leaked from the electrode into the electrolyte will reduce electrolyte. The mechanical and chemical properties of nanometer thick SEI, especially on the anode side, is critical to battery performance and durability. For Li anode, one of the unsolved challenges is Li dendrite growth. The formation of Li dendrite will consume Li metal which then leads to capacity loss and most importantly, can cause short circuit of the whole battery. Due to 7 the high activity of Li and the continuous formation of the surface layer (SEI layer), the surface morphology becomes rough, mossy, and dendritic during cycling. The morphology of Li dendrite is impacted by its complicated chemical environment and stress state. To restrain the formation of Li dendrite growth, it is important to understand the chemical-mechanical coupled effect on the morphology evolution and growth mechanism of Li metal [62,63]. Figure 1.4 Computational and experimental techniques for a variety of length and time scales [68] Novel experimental techniques were developed to measure the mechanical properties and characterize the deformation mechanisms of nanostructured materials with complex structures and coupled chemical reactions. Nano-indentation and atom force microscope (AFM) can accurately select a small region of the sample to perform indentation, pressing, bending, and twisting deformations. The measured results can be combined with the finite element method (FEM) 8 modeling to extract mechanical properties such as yield strength and modulus on nano-structured material. [64] However, due to the difficulty in conducting experiments on nano-structured material, the reported experimental measurements on their mechanical properties showed significantly scattered values [65–67]. Figure 1.4 shows different characterization techniques to observe features and investigate deformation mechanisms at different length scales. Among them, the transmission electron microscopy (TEM) and its branching techniques can achieve the atomic level spatial resolution during in-situ observation and environmental effect. The high-resolution TEM (HRTEM) observation on crystal lattice is first pioneered by Iijima [69] and Cowley [70] in the 1970s. After that, the HRTEM has been used to reveal the atomic structure of grain boundaries [71], dislocation core structures [72], and surface reconstructions [73]. Hashimoto et al. developed the environmental TEM (ETEM) which further extended the experimental window from vacuum to various gas environments[74]. After the year 2000, ETEM techniques grows rapidly that enabled its capability to observe the chemical-mechanical coupled phenomenon of nano-structured materials, such as carbon nanotubes and whiskers, in a controlled chemical environment that resembles the true experiments [75,76]. Li et al. implemented ETEM to study the degradation effect of hydrogen gas on the protective thin aluminum oxide layer. The effect of hydrogenation on metal/oxide interface detachment is shown in Figure 1.5. They proposed that the hydrogen- induced interfacial failure is due to the formation and recombination of proton interstitials and vacancies which have higher mobility at elevated temperature [77]. Another example is the impact hydrogen on dislocation movement in aluminum studied by Xie et al. [78] By observing the dislocation response under cyclic loads using ETEM, they concluded that the notorious hydrogen 9 embrittlement effect is due to the hindering effect of hydrogen-vacancy complexes rather than the previously proposed hydrogen interstitial – dislocation interaction. Figure 1.5 Comparison of the heating effect on (a,b) hydrogen-free and (c,d) hydrogenated pillars using ETEM. After annealing, the dislocations escaped from the surface for both pillars. But the hydrogenated pillars also formed a cavity between metal and surface oxides on the top of this pillar. Although many advanced experiment techniques such as TEM, HRTEM, ETEM [79] are developed to directly characterize the materials’ structure and morphology at the atomic level, only samples with the stable condition under electron beam can be examined [80–83]. These techniques do not apply to the chemically reactive and electron beam sensitive nano-structured material. With the state-of-art cryo-electron microscopy, the static atomic level information can be collected on electron beam sensitive materials such as Li nanowire since the instrument operation and sample transfer is under cryogenic temperatures [84]. Overall, due to experimental limitations, many deformation mechanisms for nano-structured materials are unknown. 10 Currently, the major challenges to understanding the deformation mechanism of complex nano- structured material are the complexity of structure and the chemical-mechanical coupled effect. Experimentally, it is difficult to isolate the contributions from individual components in a complex nanostructure, where multiple structure levels exist. While some experimental techniques are more sensitive to chemical change and others are more precise for mechanical responses, capturing chemical-mechanical coupled phenomena in nano-structures remain challenging. However, it can be envisioned that the interaction between nano-structured material and the environment is inevitable and are often enhanced because of the high surface and interface reactivity. The chemical reaction with the environment can have coupling effect with the mechanical properties and even alter the deformation mechanism. To overcome and compensate the experimental limitation, we proposed and implemented the multiscale modeling techniques to gain insight into the chemical-mechanical coupled atomic level deformation mechanism of complex nano- structured materials. As shown in Figure 1.4, the modeling methods covered the length scale of nano-structured material and have the capability of modeling mechanical-chemical response with the desired accuracy. Therefore, the focus of this thesis is to explore the nano-scale deformation mechanisms of complex nano-structured materials using predictive modeling techniques at the electronic and atomistic level. Three examples are given including the high strain rate toughening mechanism of nacre, the aluminum oxide bifilm fracture mechanism, and the lithium nanowire morphology change in the oxygen environment. While the nacre and Li-metal examples are motivated by experimental observations, the mechanical properties of aluminum oxide bifilm formed during aluminum casting have never been measured. The advantage of predictive modeling is to calculate the material properties without empirical experimental parameters. To simulate such a broad range of 11 materials and deformation mechanisms, new simulation protocols will be developed in this thesis. These modeling results not only provide new understanding of the nano-materials deformation behaviors but also give directions to further improve the materials performance, due to the predictive nature. Overall, these simulation methods open up a new forefront for modeling of nanomechanics and gaining mechanistic understanding. 1.2. Multiscale modeling techniques on deformation mechanisms Complementing traditional experimental techniques, multiscale modeling techniques including quantum mechanics methods (QM), forcefield based molecular dynamics (MD) and finite element methods (FEM) are used to understand the fundamental deformation mechanisms associated with hierarchical and complex structures. As shown in Figure 1.6, multiscale simulation techniques can be implemented to different length scale. Extremely High accuracy results of total energy can only be achieved at the lowest length scale. The accuracy is compromised as increasing of the system size since larger scale structures are approximated by a series of assumptions to reduce the complexity. The higher-level results can be calibrated using lower level information such that the important physical properties are reproduced at a satisfied accuracy. Meanwhile, new constitutive equations, material properties, physical insights provided by the lower length scale can be carried out to the higher length scales. Here, the strength and limitations of each modeling method will be introduced along with example deformation mechanisms they have revealed. 12 Figure 1.6 The multiscale simulation techniques for materials modeling at different length scale and time scale. Methods applied to the larger system usually need to compromise on calculation accuracy 1.2.1. Individual Deformation mechanism study using DFT At the lowest length scale, the quantum mechanics based first principles method is used. The total energy and electron structure of a system can be accurately calculated by solving Schrodinger’s equation (SE) without any empirical parameters. The most stable structure can also be obtained by minimizing the total energy with respect to the atomic positions. Due to the difficulty in solving the SE for systems with many-body electron interactions, a series of assumptions are introduced to simplify SE. A straightforward method is the Hartree-Fock method (HF) which ignores the electron correlation term in the Hamiltonian[85]. Many Post-HF methods are developed such as 13 configuration interaction (CI) [86], coupled cluster (CC) theory [87] and Møller–Plesset perturbation (MPn) theory [88]. These methods can be applied to small molecules and achieve highly accurate results. However, the computation cost is too high to be applied to meaningful structures in material science. Density functional theory (DFT) method is a more scalable and computationally less costly quantum mechanical method for total energy and electron structure calculation comparing with previously mentioned wavefunction-based methods. It is mainly based on two theorems from Hohenberg and Kohn [89,90]: The ground-state wave function is a unique functional of electron density and a defined energy functional is minimized by the ground-state electron density. DFT method, especially the plane-wave DFT enables the study of bulk crystals with relatively high accuracy. More detailed descriptions of theories for DFT can be found in references [91–93]. Some individual deformation mechanisms, involving hundreds to one thousand atoms, can be accurately studied using DFT method. For metals, the formation of stacking fault and twinning are important basic deformation mechanisms along with dislocation nucleation and interaction. This is often discussed in FCC metals, where the close-packed (1 1 1) planes packed in a sequence of …ABCABC…in the perfect crystal. Under the external stress during the deformation, the close-packed planes may slide along <112 > directions and lead to a “fault” in the stacking sequence such as ABC|BCA stacking. The energy difference between the faulted structure and the perfect structure is called the stacking fault energy (SFE). The total energy change as a function of sliding distance is called the generated stacking fault energy (GSFE) curve [94]. The GSFE can be evaluated by the DFT method with a slab model, and the energy barrier can be calculated using the nudge elastic band (NEB) method [95]. If a second sliding event happens on the adjacent layer to the first fault plane, it forms a stacking sequence of ABC|BCA, which is called deformation 14 twinning since the stacking sequence is mirrored with respect to the twin boundary. On both sides of the twin boundary, their stacking sequence in the perfect crystal is maintained. The energy change in this process can also be accurately calculated using DFT method [96,97]. By comparing the activation barrier of forming a stacking fault or twin, the preference of forming a stacking fault (along with a full dislocation) or a twin in an FCC metal can be predicted by an indicator, twinning tendency [98]. From this indicator, the competition between dislocation emission and deformation twins can be predicted. The GSFE also reveals other important features during deformation such as dislocation migration barrier, fracture energy, and strain energy. However, the concept of GSFE has not been widely used to analyze deformations in ionic materials. Dislocation nucleation, motion, and interaction are the most critical deformation mechanisms for plastic deformations in metallic materials. The dislocations can interact through elastic fields created by distortion of crystal lattice from a dislocation core. So the actual core structure is essential information towards prediction of the elastic fields and dislocation interactions. In f.c.c metals, a full dislocation can reduce its energy by splitting into two partial dislocations. The splitting distance is determined by the stacking fault energy. In high stacking fault energy materials such as Al, the core radius can be quite small and hard to observe in experiment [99]. So the DFT method has been used as a powerful tool to model the dislocation core structure [100–102]. The information collected from the correct prediction of dislocation core structure can be used as input of the plasticity modeling methods at larger length-scales [103]. However, due to the size limitation, so far DFT calculations are limited in identifying dislocation core structure and calculating solute-atom and dislocation interactions. Overall, the DFT method is highly accurate but computationally intensive. Currently, the system size is usually limited to hundreds to a 15 thousand atoms. Also, DFT can directly calculate only ground state energy which means the dynamics at finite temperature and stress-field are often ignored. 1.2.2. Deformation process and mechanisms studied using the traditional forcefield based molecular dynamics To predict mechanical properties at a larger length scale and finite temperature, the forcefield based molecular dynamic (MD) method need to be implemented. Due to its high scalability for parallel computing, MD simulations can simulate systems with more than one million atoms [104]. Thus direct simulation of deformation of nanostructures similar to experimental settings can be achieved [93]. Figure 1.7 Typical terms contained forcefield for organic materials. The total energy consists of bond stretching, bond angle bending, bond twisting, van der Waals interactions, and Coulomb interactions. [105] 16 In classical forcefield method, the interactions between atoms are described by a set of analytical equations and the atomic motion is governed by Newtonian physics. The analytical equations summarize the contributions of atomistic interactions explicitly based on the understanding from quantum mechanics or experimental observations. Different force fields have been tailored for different class of materials based on their bonding nature. For metallic systems, the embedded atom method (EAM) captures the shared electrons in metallic bonds [106]. For ionic systems, Bush potential [107] uses a core-shell model that consists of a core and shell carrying opposite amount of charge and connecting with a spring to describe polarization. For organic materials, many famous forcefields were developed, such as AMBER [108–110], CHARMM (Chemistry at HARvard Molecular Mechanics) [105], Universal force field (UFF) [111–114], and DREIDING forcefield [115]. The total energy in these force fields consists of bond stretching, bond angle bending, bond twisting, van der Waals interactions, and Coulomb interactions, as shown in Figure 1.7. The parameters in the forcefields are developed by fitting various material properties from experiments and energies (sometimes forces) provided by DFT calculations. So intrinsically it is not guaranteed that the models can correctly capture all material properties. Most forcefields have their applicable systems with limited transferability to other systems. For example, embedded atom method (EAM) [106] can only be applied to a metallic system because of the delocalized electron assumption. In practice, the accuracy of a forcefield heavily relies on the optimization algorithm for parameter searching, the training data to be fitted and experience from the researcher. Formulation of the forcefields and more information can be found in references [91,92]. Although the forcefield based method is less accurate in energy calculation than DFT, it has the capability to model larger systems and properties at finite temperature. Embedded atom method 17 (EAM) has been widely used to investigate the plastic deformation mechanism in metals in the form of single crystals, bicrystals, nanowires and polycrystalline nanocrystals [116–119]. Deformation mechanisms involving dynamics, such as dislocation sliding [120], dislocation cross- slip under stress field[121], grain boundary migration [122] and grain boundary sliding [123], can be carefully studied. For example, the deformation mechanism transition on metallic nanowire induced by cross-section area difference is investigated by Lao et al. using EAM potential. It was observed that the cross slip is responsible for the plastic deformation of nanowires with large cross- sectional area while the crystalline slip is dominating the deformation mechanism of small cross- sectional area nanowires [124]. Diao et al. proposed that for small cross-sectional area nanowires, the phase transformation and crystallography reorientations may also happen [125–128]. With the help of large scale EAM based MD simulation, Yamakov et al. explored the low-temperature deformation mechanism map for f.c.c metal polycrystals with different grain sizes [129]. As shown in Figure 1.8, this map correctly captures the deformation mechanism transition from dislocation- based to grain-boundary-based deformation mechanism with the decreasing grain size. It showed the complete extended dislocations or partial dislocations dominate deformation in the region I and region II. Grain boundary related mechanism dominates the deformation in region III. It provided direct evidence and explanation for the transition from Hall-Petch to inverse Hall-Petch relationship, as observed in Figure 1.1. However, the EAM potential applies only to metallic systems such as pure metals and alloys. It lacks the ability to incorporate chemical reactions into the deformation process. To study nano- structured material with complex structures and chemical-mechanical coupling effect, more advanced forcefield, such as ReaxFF reactive forcefield, needs to be implemented. 18 The classical forcefield cannot model chemical reactions because the accurate prediction of atomic charge change requires quantum mechanics approaches. In classical forcefield method, for systems with less variation on charges (e.g., ionic system), fixed formal charges are often assigned to each element and remain constant during the simulation [130]. For metallic systems, the embedded atomic method (EAM) treat electrons as uniformly distributed electron fields [106] without defining explicit atomic charge. The classical force fields for organic materials, such as AMBER, CHARMM, UFF, and DREIDING, all keep covalently bonded atom pairs bonded during simulations. These force field model cannot describe covalent and ionic bonding change during deformation and chemical reaction. Thus, it is difficult to observe the impact of environment on the deformation process with these classical forcefields. Figure 1.8 Proposed deformation mechanism map. σ/σ" is the normalized stress; r(cid:19)/(cid:10) is the inverse of normalized grain size. [129] 19 To model the chemical reactions during mechanical deformation, the bond order (typically for covalent bonds) and atomic charges need to be calculated dynamically. Different than the classical forcefields, the forcefield with the ability to model chemical reactions are called reactive forcefield. Currently, successful reactive forcefield includes ReaxFF developed by Adri van Duin[131] and COMB developed by Susan Sinnott[132]. The ReaxFF has been successfully implemented in hydrocarbon systems[133], ceramics[134], metal-and-oxide, and many other more systems[135,136]. With the development of ReaxFF reactive forcefield, it is possible to study the chemical- mechanical coupling effect on mechanical properties and deformation mechanism of materials using MD simulations. The first ReaxFF parameters that can capture both metallic and oxide was developed in 2004 by Zhang et al. for studying the non-wetting to wetting transition of an aluminum droplet on the #-Al2O3 surface [137]. Sazgar et al. simulated the fracture behavior of Al/Al2O3 interface with different termination and orientation using ReaxFF. It is demonstrated that during the deformation, the ReaxFF can correctly capture the high tensile and shear strength of O- terminated interface which is due to higher bonding strength of Al-O bond. This is in consistent with more accurate DFT calculation on fracture energy of Al/Al2O3 interface [138]. They also observed two deformation mechanisms in shear mode at different orientation: the planar sliding and amorphous transformation [139]. Sen et al investigated the oxidation effect on the tensile deformation of the aluminum nanowire using ReaxFF. The simulation showed that the formation of oxide shell can decrease the yield strength by increasing the umber of dislocation nucleation sites and lower the activation volume as demonstrated in Figure 1.9 [59]. Also, the amorphous oxide shell showed superplasticity due to the dynamic healing effect of enhanced Al/O diffusion in the O2 environment. This mechanism predicted the ductility enhancement of aluminum 20 nanowire in an oxygen environment below a critical strain rate. Later on, this superplasticity phenomenon of alumina was observed in experiments conducted by Yang et al. [60]. Not only aluminum, the mechanical properties of nickel nanowire is heavily impacted by oxidation. Aral et al. studied the deformation mechanism of nickel nanowire in O2 environment. The oxidation layer have similar dislocation nucleation effect which leads to early initiation of plastic deformation and lower yield stress on oxide-coated Ni NW with a smaller diameter (~5.0 nm) [139]. Figure 1.9 Stacking fault formation in deformed aluminum nanowire in (a, b) vacuum and (c, d) O2 environment. Color bar indicates the central symmetry parameters. More dislocation nucleation events are observed in the O2 environment. The impact of lithiation reaction on mechanical properties of silicon has also been studied using ReaxFF. Silicon is one of the most promising anode material candidates for Li-ion battery due to its high theoretical specific capacity (3750 mAh/g, ten times higher than graphite). However, it application is limited by the 300% volume expansion during lithiation process which may lead to 21 dramatic structural change and fracture of electrode and SEI layer. The brittle to ductile transition of Si during the lithiation process is studied by Wang et al. using ReaxFF. The LixSi slab showed brittle fracture when x is less than 2.5. When x is larger than 2.5, the fracture tip is blunt because Li can form chemical bond with Si during the fracture process. To mitigate the mechanical failure issue, many techniques have been developed such as applying a protective coating layer on Si [140–142]. By implementing ReaxFF, Kim et al. simulated the lithiation process of Si-core/Al2O3- shell and Si-core/SiO2-shell nano-structured cluster in the anode of Li-ion battery [143]. It showed that Al2O3 coating could provide more resistance to the volume expansion of the Si core, and avoid stress concentration which potentially decreases the capacity loss. This stress gradient difference is due to the different chemical reaction nature between lithiating Al2O3 and SiO2. This thesis will use ReaxFF to simulate the formation and fracture of oxide bifilms in cast aluminum and the Li growth morphology in different oxygen partial pressure. Both examples require designing new MD protocols to mimic the experimental conditions. 1.2.3. Deformation mechanism studied using FEM At the macroscopic level, the material’s mechanical responses to various load conditions can be predicted from finite element methods. In the FEM method, the material is subdivided into small domains (finite element), whose deformation is continuous and is governed by constitutional partial differential equations. With known boundary value conditions applied to each domain, the numerical solutions of internal stress can be solved for the deformation process. FEM can also resolve the deformation details of microstructures, with the information of crystallographic orientation, elastic modulus, yield strength, and plasticity. However, the FEM method relies on the known constitutive relationship, which must be provided as input information. It can model 22 macroscopic deformation behavior but lack of the ability to discover new deformation mechanism at nanometer length scale [144,145]. This thesis will not utilize FEM simulations. However, one aim of the predictive modeling is to provide material properties as inputs to simulate the macroscopic performance of materials at components level using FEM methods. This connection will be constructed for the fracture energy and strength of aluminum-oxide bifilms in virtual casting simulations at the FEM level. 1.3. Motivation and Dissertation Outline This thesis will design and apply predictive modeling techniques, including DFT method and forcefield based MD methods, to understand deformation mechanisms of complex nano-structured material and describe the chemical-mechanical coupled interactions. The thesis is organized as follows: In CHAPTER 2, atomic level deformation mechanism under dynamic loading condition in nacre is proposed. This study was motivated by recent experimental work by Huang et al., who reported higher strength and the same level of toughness of nacre under dynamic tensile loading conditions at a strain rate about 103 s-1 [51]. In the specimens fractured under high strain rate, they observed deformation twinning or stacking faults with HRTEM [146]. These observations raised a series of interesting questions: what is the toughening mechanism under dynamic loading condition at the atomic level? Why stacking fault and deformation twins can be formed in biomaterial nacre but not the mineral aragonite CaCO3, although they share the same crystals at the nanoscale. To address these questions, we implemented DFT method combined with the nudged elastic band method (NEB) to calculate the generalized stacking fault energy (GSFE) curve and twinning tendency in aragonite CaCO3. Although GSFE has been routinely computed for FCC 23 metals, it has not been used to analyze twinning and stacking fault formation in crystal structures as complex as CaCO3. Therefore, we carefully designed the slip planes and twin formation according to HRTEM reported crystalline orientations. We also applied different mechanical constraints normal to the sliding direction to mimic the condition in a bulk crystal or nano-grains surrounded by soft bio-polymers. It was discovered that a critical difference between biomineral and mineral material is the strain relaxation during sliding, which can change the competition between twinning vs. fracture. Another observation was the highest ratio of the unstable stacking fault energy and the stacking fault energy in aragonite compared with a broad range of metals, semiconductors, and oxides. We further explained the underlying physics for this high ratio is the multi-neighbor shared ionic bonds and the unique relaxation process during sliding in the aragonite structure. This nanometer level deformation mechanism triggered by the high strength under high strain rate works together with other energy dissipation mechanisms to achieve the overall extraordinary mechanical performance of nacre. Overall, the understanding of the orchestrated multiscale deformation mechanisms in nacre brings new insights for the design of high toughness materials at high strain rate. In CHAPTER 3, the formation and deformation mechanism of oxide bifilms in cast aluminum is investigated. This study was motivated by a long-lasting question in the aluminum casting process. During the aluminum casting process, the liquid aluminum will be exposed to air inevitably and then form aluminum oxides rapidly [61]. Next, the aluminum oxides will be folded by the turbulence of injected liquid aluminum and trapped into the final casting product. This folded alumina thin film is called oxide bifilm in the casting industry. The initially formed bifilm is called “young oxide” with an amorphous structure. After holding in furnace for serval hours, the “young oxide” will undergo a series of phase transition from amorphous to gamma phase and 24 eventually become “old oxide” with alpha phase. Furthermore, the surface of bifilm can react with CO2 and H2 from the surrounding environment or dissolved in the aluminum matrix. So the actual bifilm surface may be covered by contaminations such as aluminum hydroxide or aluminum carbonate. The thickness of oxide bifilm is only a few nanometers [57,58], but it is highly detrimental to the casting product’s mechanical properties such as cycle life and toughness [147]. Due to the difficulty to characterize this ultrathin bifilm, complex crystallography and its complex chemical reaction with the environment, the fracturing and deformation mechanism is still unknown. Meanwhile, at the macroscopic level, virtual casting based on the FEM method has been implemented to simulate the distribution of bifilms in the cast aluminum components [148,149]. The result agrees well with x-ray radiography experimental techniques. By implementing immersed element-free Galerkin method, the resistance of bifilm to deformation can also be modeled. This model can be used to investigate the oxide bifilm folding process during solidification [150]. However, to simulate the impact of oxide bifilms on the mechanical property of the final cast components require the knowledge of the mechanical strength or fracture energy of the oxide bifilms. To fill this gap, we designed a protocol to form oxide bifilms starting from modeling the oxidation process of aluminum. We built the oxidized amorphous layers to mimic the “young” oxide and α- Al2O3 layer to represent the “old” oxide. We further included the hydrogen termination in “old” oxides. The folded oxide bifilms were then subjected to uniaxial tensile MD simulations. The oxidation and deformation processes were simulated by ReaxFF reactive forcefield based MD. Overall, we explored the impact of surface contamination, crystallography evolution and folding process on the mechanical properties of the ultra-thin Al oxide bifilm trapped in the aluminum 25 casting products. The results showed that the toughness and fracture strength of oxide bifilm is sensitive to the formation of hydroxyl group surface contaminations. We also proposed a bridging method to transfer the MD predicted fracture energy and strength to FEM simulations. The understanding of deformation and fracture mechanisms of bifilm under the influence of chemical and structural evolution provided insights into a more rational design of casting procedure to achieve high-quality casting product. In CHAPTER 4, the morphology change Li grown in different oxygen environment was simulated, and the competition between oxidation rate and Li growth rate was discussed. This simulation was motivated by recent experiments by Alexander Yulaev et al. To examine the Li dendrite formation process closely, a carbon/LiPON/LiCoO2 solid-state battery was synthesized, and the Li-growth under a voltage bias on the back of the carbon anode was observed by in-situ SEM [151]. A bias is applied to the solid-state battery such that the Li moves from LiCoO2 layer to the carbon layer and grows out from the back of the carbon layer (the top surface from view of SEM). The result showed that the Li nanowire growth morphology is heavily influenced by oxygen partial pressure and the electric current density. Under the ultrahigh vacuum condition, the Li growth on anode showed 3D mode. As oxygen partial pressure increased, the growth mode become 1D nanowire growth at 5.7×10(cid:22)( Pa. However, if the partial pressure increases to 4.8×10(cid:22)+ Pa, the growth mode become 3D again. The growth mode of Li nanowire is closely related to the notorious Li plating and dendrite formation problem on the anode that damages battery performances and may cause safety issue [152]. So, the understanding of Li growth mechanism under the impact of the oxidation reaction and electric current density coupling effect is critical to better and safer battery design. It is proposed that the morphology transition is due to the difference of Li nanowire structure under different oxygen partial pressure. 26 To clarify the competition between Li-growth rate and the oxidation rate, MD simulation based on ReaxFF reactive forcefield was performed. The atomic model showed the competition between oxidation rate and Li nanowire growth rate (controlled by the current density) changes the microstructure of Li nanowire which leads to different growth mode. From the MD observation at the atomic level, we proposed the formation mechanism of Li nanowire with different oxygen partial pressure and electric current density which explains the 3D-1D-3D morphology change of Li nanowire formed on top surface of the carbon layer in the stacking battery. When the oxidation rate is in deficient comparing to the growth rate of Li nanowire, the formed oxide layer is too thin to support the Li-1D growth. When extra oxygen is present, the oxide layer is too thick, so the 1D growth is suppressed. This result demonstrates the importance of trace amounts of preexisting or ambient oxidizing species on lithiation processes in solid-state batteries. Understanding the impact of chemical reaction on Li nanowire deformation mechanism and morphology evolution also provided us insights on the formation of the SEI layer in Li-ion battery. In CHAPTER 5, we discuss why ReaxFF might over-ductile ionic materials. Despite the successful simulation of chemical-mechanical properties with ReaxFF, there might be an intrinsic limitation of the current ReaxFF, which dynamically determines the charge with the qEq method. The qEq charge transfer method is known to have an unphysical long-range interaction and residual charge when atoms are far apart. This will bond two atoms with Coulombian interactions even they are far away to each other. Thus the bond breaking may not be appropriately described, which leads to ductile fracture even in brittle materials. To mitigate the charge transfer problem and have an accurate description of the charge transfer process during large deformation and fracture process, we applied the ACKS2 charge transfer method to ReaxFF reactive forcefield. The ACKS2 method was one of the recent techniques to eliminate the residual charge 27 problem known for qEq [153], and it has been implemented in PuReMD code by our collaborators in Dr. Aktulga’s group [154,155]. To compare the impact of the charge transfer process on the deformation process between qEq and ACKS2 method, a Li2O slab model with a sharpy test notch was built. This slab was subjected to a tensile test using both qEq and ACKS2 method. The ACKS2 parameters are fitted with high-level quantum chemistry calculations. The results showed that qEq method predicted ductile fracture in Li2O. The tensile test is simulated using ACKS2 method with exactly the same procedures. 28 CHAPTER 2. Chemical bonding origin of deformation twinning formation in the biomineral aragonite (Part of this chapter is adopted from Ref. [156]) 2.1. Summary Deformation twinning rarely occurs in mineral materials which typically show a brittle fracture. Surprisingly, it has recently been observed in the biomineral aragonite phase in nacre under high rate impact loading. In this chapter, the twinning tendency and the competition between fracture and deformation twinning were revealed by first principles DFT calculations. The ratio of the unstable stacking fault energy and the stacking fault energy in orthorhombic aragonite is hitherto the highest in a broad range of metallic and oxide materials. The underlining physics for this high ratio is the multi-neighbor shared ionic bonds and the unique relaxation process during sliding in the aragonite structure. Overall, the unique deformation twining along with other highly coordinated deformation mechanisms synergistically work in the hierarchical structure of nacre, leading to remarkable strengthening and toughening of nacre upon dynamic loading, and thus protecting the mother-of-pearl from predatory attacks. 2.2. Introduction 2.2.1. Hierarchy structure of nacre Nacre, also known as mother-of-pearl, is the shiny inner layer of the mollusk shells. The main component of nacre is aragonite (a metastable polymorph of CaCO3 with PMNA symmetry), 29 which occupies 95% of the total volume. The other 5% (volume) is biopolymer interlayers between the aragonite platelets and the aragonite nanoparticles inside the platelets. In the effort to design new materials, Mother Nature is always a good guide for material scientists to study and mimic. The shells of mollusks, such as nacre (mother of pearl), are excellent examples of toughening design for us to learn from[156]. Nacre is the shiny inner layer of the mollusk shells. 95 vol% of nacre is composed of aragonite, which is a meta-stable phase of CaCO3 at room temperature. The other 5 vol% is a biopolymer[157]. Nacre has a hierarchical structure. As shown in Figure 2.1, the first order structure is 1μm thick layered structure. The layered structure is divided into small platelets with a diameter of 1-5 μm and thickness of 200-500 nm, which is the second order structure. Between platelets, a biopolymer is filled in. Furthermore, the second-order platelets consist of pseudo-hexagonal platelets that are 30-80 nm width [19,146,158]. Previously, the third order platelets are treated as single crystal materials that have the same mechanical properties with bulk calcium carbonate. With the development of TEM and SEM technology, later studies showed that they consist of smaller structures [19]. 30 Figure 2.1 Top: the hierarchical structure of nacre[52]. Bottom: The first order structure is a layered structure. The second order structure is platelet that’s further divided into smaller sub- plates [159]. The “Brick and Mortar” structure in nacre shell can be observed by microscopy directly. As shown in Figure 2.2 [160], multiple layers of aragonite layers are piled up and interlocked with each other. Between each platelet, there is a 20-30 nm gap that filed with β- chitin, a bio-polymers, acting as the “glue” to connect each platelet. Zhi-Hui Xu et al. measured the elastic modulus of the bio- polymer to be 11±3 GPa using AFM tip and finite element method[64]. What is more, Schäffer et al. observed nanopores on the biopolymer matrix indicating the platelet are also connected by mineral bridges [161]. F. Song et al. [162] observed that on each of the individual platelet, small pillars with 45 nm diameter are connecting adjacent platelets which validated the existence of mineral bridge. A schematic illustration of the platelet structure is shown in Figure 2.1 where the polygonal platelets are the “brick” and bio-polymers between the gaps are “mortar.” 31 Figure 2.2 The TEM imagine of platelet structure in nacre shell [160] A further examination of nacre shell showed that on the top of each pile, we could see several nanoparticles as seeds which also indicating the growth process of each platelet [160]. Gangsheng Zhang et al. directly observed that at first, aragonite nanoparticle self-assembled into a core-shell structure and then form a dome-shaped platelet. Previously, the platelets are treated as single crystal that has the same mechanical properties as bulk aragonite because the small angle diffraction (SAD) measurements of platelet clear showed a single crystal pattern [163]. However, the self-assemble growth mechanism implies that each platelet consists of many small grains instead of being a single crystal. 32 Figure 2.3 A schematic illustration of nacre shell structure including the platelet, biopolymer and mineral bridges [162] 2.2.2. Toughening mechanisms of nacre under quasi-static loading Compared with the mineral aragonite, nacre’s toughness is a thousand times higher, and its mechanical strength is two times higher[163–165]. Such attractive physical properties have been intensively investigated and ascribed to nacre’s hierarchical structure, where the brick-mortar microarchitecture can regulate crack propagation in a zigzag manner[51,156,166,167]. Yet, little attention has been paid to the protection mechanism activated in nacre to maintain its mechanical integrity in the event of penetrating impact such as predatory attack, tidal current[168] and so on, which often occurs at high speeds. 33 Figure 2.4 Illustration of feature of platelet protruding in interlock mechanism [53] The platelets are locked at the circles to prevent relative displacement and rotation Under static loading, the generally accepted toughening mechanisms are the following. The inorganic platelets are the most fundamental component of nacre shell. It undertakes most of the load and absorbs large percent of the energy during inelastic deformation. The hierarchical platelets directly improve the toughness by guiding the propagation of crack and allowing slight inter-platelet slip to increase the ductility. Other than these, there are still many other features the shell evolved to increase toughness. The interlock mechanism proposed by Katti et al. [53] is one of the most representative mechanisms based on the platelet structure. The edge of platelet protrudes out about 20-30 nm and locks the rotation and pull out of next platelet (see Figure 2.4) when the tensile load is applied parallel to the surface of shell. A FEM model is established in order to evaluate the impact of interlock. The results showed the interlock mechanism can improve 34 the ultimate tensile stress from 14 MPa to 50 MPa [169]. But this is still not comparable to a real shell which has a 170 MPa tensile stress[156]. Figure 2.5 The short/long molecule model [54] (a) Schematic drawing of the biopolymer chain between two platelets. The long molecule chains on the top can be folded to modules and become short chain on the bottom. (b) The response of a short/long molecule chain under extension. The elongation can be increased while maintaining the flow stress by opening up modes on the molecule chain. The biology effect of biopolymer is to regulate the deposition and assemble of inorganic platelet at every level. But in fact nacre shell contains far more bio-polymer than the need for regulating platelet growth. This indicates that although the bio-polymer only occupied about 5% volume fraction, it played an important role in the mechanical properties [170]. The main component of biopolymer is chitin and there are also a small fraction of proteins [157]. The stickiness of biopolymer is the first property people focused on. In the theory of composite materials, adding a polymer matrix will greatly improve toughness. This is also verified by synthesizing model 35 materials. Furthermore, Smith et al. [54] stretched the exposed bio-polymer on a freshly cleaved nacre and discovered the step-wise manner in the elongation process. They proposed a short/long molecule model to explain the step-wise behavior (See Figure 2.5). The short molecule can sustain high strength when the strain is low while after the break of short molecules, the long molecule can still absorb energies. Also, instead of breaking all the molecule at the same time, short ones are sacrificed first to dissipate the stored energy into heat. During the quasi-static deformation, a large amount of crack goes through the bio-polymer instead of the inorganic platelet because of its relatively low modulus [51]. Most energies are absorbed by bio-polymers in quasi-static deformation. The crack deflection mechanism proposed that the crack tend to propagate along different directions in adjacent layers which will prolong the crack path, absorb more energy and generate a tortuous fracture surface. As the deformation keep going, platelets near the fracture surface will be pulled out eventually. In the pullout process, bio-polymer matrix bridging between two platelets is observed by R.Z. Wang et al. [171]. Implementing finite element method, Katti et al. [172] estimated the elastic modulus of the interface between biopolymer and inorganic platelet to be 15 GPa while Xiaodong Li et al. [64] measured the overall elastic modulus of biopolymer matrix is 11 GPa. This ensured the binding in the interface during stretching of biopolymer and indicates more energy can be absorbed during the pullout process. 2.2.3. Experimental observations under dynamic loading Recent uniaxial compression tests showed that the fracture strength of nacre under the impact loading (strain rate: ~103/s) is approximately three times higher[146,173] than that under quasi- static loading (strain rate: 10-3/s). The actual strength was increased from 170 MPa to 500 MPa as shown in Figure 2.6 (a). At the same time, the compressive plastic strain remains the same or even 36 further improved[51]. Moreover, its toughness (characterized by the strain energy on the true stress-strain curve) increased from 2.4 to 6.1 mJ/m3 as the loading condition changes from quasi- static to dynamic. Figure 2.6 Experimental observation after deformation of nacre (a) Stress curve of nacre shell under different deformation rate and the HRTEM imagine of specimens after deformation. After quasi-static deformation, only (b) grain rotation is observed. But after dynamic deformation, (c) twinning and (d) partial dislocation emission can be observed. [146] High-resolution TEM images are collected from the samples after deformation at different conditions. After quasi-static deformation, twinning and partial dislocation are observed in nacre as shown in Figure 2.6 (c) and (d). Under quasi-static loading condition, only grain rotation is observed in HRTEM as shown in Figure 2.6 (b). No deformation twin or partial dislocation was 37 found. The grain size is only about 20 nm in diameter which is only observable using HTREM. This is the smallest sub-structure of nacre shell and different than previous understanding that the platelet is a single crystal structure. These are identified as deformation twins rather than growth twins. Previous literature has extensively documented that growth twins are ubiquitous in cross lamellar seashells[174,175] but seldom in gastropod class. Typically, the length scale of growth twin crosses throughout the tablet [176] and deformation twins are confined within an individual aragonite nano-grains owing to the extremely small activation volume. This is contrary to the nano-grain deformation and rotation that occurs under quasi-static loading[49,177]. Typical engineering materials, such as metals and polymers, become stronger but more brittle under impact loading, since many deformation mechanisms, i.e., chain relaxation, dislocation climbing, grain boundary sliding et al., do not have enough time to be activated at high strain rates [178–180]. Only a few composite materials have been engineered to enhance both the strength and ductility under dynamic loadings [178,181]. The record amount of 2.5-fold increase in toughness and 3-fold amplification in strength holds a fascination to guide the design of stronger-and-tougher materials, exceeding typical engineered materials to date as shown in Figure 2.7 (a) [146,178,182]. However, the deformation mechanism and the reason for nacre’s extraordinary mechanical condition remains unknown. We proposed that the remarkable strain rate induced strengthening-and-toughening manifestation is a result of the highly coordinated deformation mechanisms synergistically working in the hierarchical structure. Under high strain rate compression, the typical polymer chain relaxation between aragonite platelets and nanoparticles do not have time to participate in deformation. Thus, 38 strengthening and toughening are primarily carried by the inter-nanoparticle deformation/rotation, in combination with the emission of partial dislocations and the onset of deformation twinning within individual aragonite nanoparticles. These nanoscale deformation mechanisms have been held responsible for the extremely small activation volume (~ a few b3). Figure 2.7 Nacre’s extraordinary strengthening and toughening and its multi-functional deformation mechanisms. (a) The normalized toughness performance of different materials under high strain rate (~103) and low strain rate (~10-3). Nacre shows dramatic toughness increase over metals, polymers and even traditional composites under high deformation rate. (b) The deformation mechanisms that can be activated at the atomic level. When the deformation rate is high, only the mechanisms illustrated in grey platelets can be activated. However, the atomistic origin of deformation twinning in aragonite under impact loading remains completely unknown. Moreover, why mineral aragonite (monocrystal form) behaves brittle, while 39 the nanoscale deformation twining/partial dislocation can be triggered in biomineral aragonite, need to be further elucidated. But due to the nano-scale length of the deformation mechanism, it’s extremely difficult to conduct in-situ experiments. In this study, we implemented atomistic modeling techniques based in density functional theory (DFT) method to present the atomistic insights into the unique twinning formation mechanisms in the biomineral aragonite phase and uncover the competition mechanism between nanoscale plastic deformation and brittle fracture. 2.3. Model Construction and Computational Methods 2.3.1. Compute the GSFE and Twinning Tendency GSFE was first introduced by Vitek et al.[183] to quantify the energy change during stacking fault formation due to atoms sliding after a partial dislocation emission in f.c.c metals. It is easy to illustrate twinning and stacking fault in face-centered cubic (fcc) crystals. Its stacking sequence of (1 1 1) plane is …ABCABC… in a perfect crystal. During deformation, if some (1 1 1) layers slide along the <112> direction, the layers shift positions accordingly, (A→B, B→C, C→A), and the stacking sequence will become …ABC|BCA…, which forms a stacking fault, or partial dislocation. If the structure keeps on shifting on the same slip plane, the stacking sequence will be restored and left a full dislocation behind. But if the structure slide on the plane next to the stacking fault plane, the stacking sequence will become …AB|C|BAC… and the sequence is mirrored with the twinning boundary, the C plane in the middle. This is called a twin structure with the thickness of 2 layers. From the fault structure, either a full dislocation or a small unit of twin structure will form depending on whether the next slide is on the same slip plane or next to the stacking fault plane. The fault energy γ per unit area can be calculated as a function of sliding displacement, u, and the fault energy profile of γ(u) is defined as GSFE landscape. The valleys on the GSFE correspond to 40 the stacking fault energy, SFE, and twin formation energy, TFE. The peaks of the GSFE represent the energy barriers to form a stacking fault or a twin, thus defined as the unstable stacking fault energy, USFE, and the unstable twin formation energy, UTFE, respectively. Initially, embedded atom potential (EAM) has been used to evaluate the GSFE curve, but the accuracy varies with the quality of the fitting of forcefield. Swygenhoven[94] calculated GSFE curve in nanocrystalline Ni, Cu and Al systems with different EAM potentials. Although the stacking fault and full dislocation formation were captured in the simulation, the SFE, USFE, and TFE calculated from Cleri–Rosato’s potential was three times higher than Mishin’s potential. G. Lu[99] compared the GSFE curve of Al calculated by density functional theory (DFT) with Local Density Approximation (LDA) with that from Ercolessi-Adams’ EAM potential. The results showed although EAM gave the general trend of GSFE, it failed to calculate the exact energy values and predict dislocation core structure. To overcome the large variation of GSFE computed by empirical potential, Siegel[184] and Wu[185] calculated the GSFE using DFT calculations and twinning tendency of Ni and Ni-based alloys. Their results confirmed the accuracy of DFT calculated GSFE and indicated alloying with transition metals will decrease SFE and increase the twinning tendency which agreed well with experimental results. In this study, the GSFE curve was calculated using climbing image nudge elastic band method (NEB)[186–189] method based on DFT calculation. NEB is a widely used method to find the transition state between two different configurations. It maximizes the force along reaction direction that connects the starting and end configuration while minimizing the force perpendicular to the reaction direction. 41 2.3.2. A short introduction to DFT method Schrödinger equation (SE) is the governing equation of classical quantum physics. The solution of the Schrödinger equation, also known as wavefunction, contains all information of the particles in a given state. But the analytical solution of Schrödinger equation can only be found in a few simple systems such as one electron in an infinite potential well and hydrogen-like systems which don’t contain the electron-electron interactions. Even for the numerical solution, the computational cost is too high to be applied to practical systems. To get an approximated result under proper time complexity, a series of assumptions are introduced to simplify the Schrödinger equation. Since the mass of nuclei is much higher than electrons, it is assumed the nuclei do not move the position and electrons respond more rapidly to the external potential in Born-Oppenheimer approximation (BO). So the wavefunction of nuclei and electrons can be separated in SE. BO leads to the time-independent, non-relativistic Schrödinger equation as shown in Equation 2.1: 1 ,(cid:28)ℏ(cid:24)2(cid:15)./0(cid:24) 02(cid:23) 1 +.3(40) 02(cid:23) 1 +..(cid:2)(40,46) 02(cid:23) 670 8Ψ=iℏ∂Ψ∂t Equation 2.1 where m is the electron mass; Ψ is the wavefunction; N is the number of particles. The three terms in the bracket are called Hamiltonian which contains kinetic energy of electrons, the interaction between electron and nuclei, and the interaction between electrons. The third term, electron- electron interaction is the most challenging part to solve since we cannot apply the separation of variables tricks in the Coulomb term. 42 To further simplify the SE, Hartree method is proposed that treat the electrons as non-interactive particles so the total wavefunction can be replaced as the multiplication of wavefunction of each electron. But this treatment ignored the anti-symmetry property of wavefunction. Fock implemented the anti-symmetry constrain by introducing the Slater determinate. In Slater determinate, swap of two particles’ wavefunction will naturally lead to a negative sign in the total wavefunction which satisfies the anti-symmetry constrain. With the introducing of Slater determinate, by applying variations theorem, we can get the Hartree-Fock equation (HF). In HF, the Coulomb term is approximated by the central field approximation which assumes each electron only interact with an electron density distribution rather than all other electrons. HF method considered the electron exchange contribution but ignored the electron correlation contribution. The exchange energy is the direct result of introducing Slater determinate. The missing of correlation energy is because the electron density is not accurate near the location of each electron. To solve for the correlation energy, many methods are developed such as perturbation theory (MP2). Currently, the most widely used method is the density functional theory (DFT). In 1964, Kohn and Hohenberg proposed the two theorems as the foundation of DFT [90]: Theorem 1: The ground state energy of a many-electron system are uniquely determined by the electron density distribution in the space Theorem 2: The electron density that minimizes the energy of overall functional is the correct electron density. The first theorem reduced the degree of freedom of SE from 3N to 3 dimensions which converted the intractable many-body problem to a tractable problem. The second theorem provided a method 43 to calculate the ground state electron density and the corresponding energy. In fact, it is a result of applying variational theory. Later, the HK theorem was further developed by Kohn and Lu Jeu Sham. They proposed the Kohn – Sham (KS) equation by replacing the interacting, whole system functional with non-interacting, separate functional [89]. All errors are collected by an extra exchange-correlation term as shown in Equation 2.2. =(cid:28)ℏ(cid:24)2(cid:15)/(cid:24)+3(>)+?(cid:24)@ A(>B) |>(cid:28)>B|+D(cid:5)EF(>) DA(>) GΨ=iℏ∂Ψ∂t Equation 2.2 In the Hamiltonian of the KS equation, the first and second term is still the kinetic energy of electrons and the electron-nuclei interaction term. The third term is also called Hartree potential, which describes the electron-electron Coulomb interaction. The last term is the exchange- correlation (XC) potential which does not have an exact formula for now. To find an approximated XC term, many approximations methods are proposed such as local density approximation (LDA) and generalized gradient approximation (GGA) [92]. LDA assumes the XC term is only related to the EXC rather than the gradient of EXC. This simple approximation leads to surprisingly well results. GGA pushed one more step by taking the gradient into consideration. In actual DFT calculation, pseudopotentials are implemented instead of all electron density distribution. It is unnecessary to calculate the exact electron density of all electrons since only shell electrons will participate in the bonding process. Also, the core electron density experiences rapid 44 fluctuations due to the nodes in wave function. This breaks down the assumption of LDA and makes it somewhat unreliable in the core region unless a large number of basis set is used. To serve our purposes, we can neglect the core electrons and only solve for the electronic structure of valence electrons. The fluctuation of core electron is replaced by a smoothly fitted curve called pseudopotential. To process the wavefunction in a computer, they are usually projected into basis sets as shown in Equation 2.3: ΨH=ΣJcJφJ Equation 2.3 where cJ is the coefficients and φJ is the basis sets. There are mainly two category for basis sets: the atomic orbitals and plane waves. For plane wave basis sets, it assumes that the nuclei are periodic perturbations to the free electrons. More accurate results are usually achieved by including a greater number of basis, but in turn, it takes more time to compute. When implementing DFT method in the material system, the atoms are usually in the periodic crystal structure. This means all atomic positions will repeat themselves after a certain distance. According to Bloch’s theorem, if the atomic structure is periodic, the solution of the KS equation must satisfy Equation 2.4: ΨH(M)=exp(iM∙4)uS(4) Equation 2.4 45 where TU(4) is a periodic function in real space. With this form, it indicates that we can solve the KS equation for each k individually. In this case, it is more convenient to solve the DFT problems in the space vector k lives in, which is called reciprocal space or k space. In practice, it is unnecessary to compute every point in k space. We can sample some important points in the k space with assigned weights and numerically calculate the integral over the whole space. Also, due to the symmetry of k space, we can reduce our sampling to the first Brillouin zone to reduce computation time. If too few numbers of k points are used in DFT calculation, the result will have a larger systematic error. So, for every DFT calculation, the number of k points must be carefully tested to ensure the convergence of final results. If we compare Equation 2.3 and Equation 2.4, the TU(>) actually contains a special set of plane waves as shown in Equation 2.5: TU(4)=.VWXYexp (Z((cid:21)+Y)∙4) Y where G is defined by unit cell vectors in k space: Y=(cid:15)(cid:23)[(cid:23)+(cid:15)(cid:24)[(cid:24)+(cid:15)\[\ with m as integers. Since it is unrealistic to calculate infinite terms over all possible G, we also need to truncate the summation to some extent. The G is related to kinetic energy in solving KS equation: Equation 2.5 (cid:5)= ℎ(cid:24)2(cid:15)|(cid:21)+Y|(cid:24) We usually refer the amount of G used in DFT calculation as cutoff energy. The choice of reasonable cutoff energy is selected by evaluating the energy of a unit cell and choosing the lowest 46 cutoff energy that satisfies the accuracy requirement. This process along with choosing of k points are part of the convergence test. In this chapter, to address the unique deformation twinning in bio-mineral aragonite, the generalized-stacking-fault energy (GSFE) curve along the sliding direction was computed from density functional theory (DFT) method implemented in the Vienna ab initio simulation package (VASP)[183,190–197]. Generalized gradient approximations (GGA) with Perdew-Burke- Ernzerhof (PBE) parameterized pseudopotential was used to evaluate the exchange-correlation function. The valance electron configurations were: C-2s22p2, O-2s22p4 and Ca-3s33p64s1. The total energy was converged to 1 meV/atom with 550 eV cutoff energy and 4×4×2 k points for the unit cell structure. The total energy of all the slab models was calculated by relaxing only atom positions while fixing the cell dimension. The energy minimization was performed with the conjugate gradient algorithm by minimizing Hellman-Feynman forces to less than 0.02 eV/Å. To calculate the unstable stacking fault energy and unstable twinning energy, the climbing image nudge elastic band method (NEB) is used in the transition state search[186–189,198]. By maximize the Hellman-Feynman forces along the sliding direction and minimize in every other direction, the saddle point is reached with a deviation of less than 0.02 eV/cell. 2.3.3. Construction of nacre slab and computation details The atomistic model for aragonite CaCO3 was constructed based on the slip systems observed in samples harvested after quasi-static and dynamic compression tests of natural nacre under high- resolution transmission electron microscope (HRTEM). The formation of stacking fault structures and deformation twinning (Figure 2.8 a) was identified and was believed to dominate the deformation process upon dynamic compression[146]. The atomic-scale lattice image (Figure 2.8b) 47 reveals that the twinning structure has a twinning boundary of (110) planes with an interlayer distance of 4.2 Å, the angle between (110) and (010) planes is 31.9° [146]. Based on the HRTEM findings, a slab model of aragonite with (110) planes was constructed. The projections of the 8- layer slab with three layers of the twinning structure are shown in Fig 2c. The dashed lines indicate the corresponding twinning planes in Figure 2.8b. The slab model is periodic along [1 10] and [001 ] directions, with a 10 Å vacuum layer added between (110) surfaces. The simulation cell dimension is 5.81Å×9.47Å×47.62Å along the x, y, and z directions, respectively. Figure 2.8 The atomic structure of deformation twinning (a) Upon impact along nacre’s cross section, deformation twinning/stacking fault takes over the deformation process, rather than the full-dislocation slipping[146]. (b) The atomic-scale measurement of white box area in (a) clearly demonstrates the twinning relationship between (110) and (11 0) planes. (c) The atomic 8-layer slab model is built according to (b) and is shown in two projections 48 2.4. Results and Discussion 2.4.1. Surface Energy of Aragonite CaCO3 The slabs with perfectly stacked (110) planes at different thickness were constructed in the middle of simulation cell so at least 10 Å of vacuum was added above and below the slab to prevent the interaction between top and bottom surfaces. The surface energy of the (1 1 0) plane is given by Equation 2.6 (cid:5)‘=a∗(cid:5)cdeU(cid:28)(cid:5)‘efg 2∗h Equation 2.6 where Es is the surface energy, N is the number of layers in the structure, EBulk is the relaxed total energy per layer obtained from the bulk structure, ESlab is the relaxed total energy of the slab model with surfaces and A is the cross section area of the simulation cell. Assuming Ebulk and ES are constant numbers, Equation 2.6 can be transformed into Equation 2.7 [199]: (cid:5)‘efg=a∗(cid:5)cdeU(cid:28)2 ∗(cid:5)‘∗h Equation 2.7 By fitting the total energy of the slab, ESlab, as a function of the number of layers, N, in the slab model, the bulk energy for one layer, EBulk, can be obtained from the slope. Plug in this value into Eq.1, the surface energy for slab with different number of layers can be calculated and the results are plotted in Figure 2.9. 49 Figure 2.9 The surface energies of different layers and their relative errors comparing with 6- layer slab model. Figure 2.9 shows how surface energy converges with slab thickness. The surface energy for 1- layer slab is much lower than the other slabs. This is due to the weak interaction between the top and bottom surfaces in the 1-layer slab. As the number of layers increases, the surface energies are converged to 491.52 mJ/m2 and for the slab with more than three layers, the relative errors comparing with 6-layers slab is less than 0.2%. This means the surface energy converges quickly with the slab thickness, thus the center of the slab is mimicking the bulk behavior and the self- interaction between the top and bottom layers is negligible. As the surface energy is converged with three layers, in order to eliminate the surface impact on the stacking fault structure, it is 50 reasonable to have three layers of atoms between the stacking fault plane and the surface. In this study, a slab model with six layers is built in the calculation of stacking fault energy. 2.4.2. Generalized SFE and Fracture in CaCO3 In the original GSFE model, sliding proceeds by single atomic layer in metals. However, as shown in Figure 2.8c, the aragonite crystal structure is much more complicated than metals. In aragonite, each (110) plane includes four CO3 2- groups and four bonded Ca2+ ions (green). The four CO3 2- groups are parallel to (001 ) plane but at different positions along the [001 ] direction, as indicated each (110) layer is equally shifted with respect to the adjacent layer by 2.65 Å along [ 1 10] by their dark blue to light blue color. In the slab with perfect stacking sequence, the position of direction. Thus, the stacking of the (110) planes in aragonite can be considered as …GABCDEFG…. The structure repeats itself every seven layers. In order to create the first stacking fault, the top four layers are shifted along the [ 1 10] direction at a distance of 4.161Å. The stacking fault structure has a stacking sequence of …BC|DEDE|FG… By shifting the top 3 layers once more, the twinning structure in Fig 2c has the stacking sequence of …DE|FEDEF|G... The calculated GSFE obtained from the 6- and 8-layer slabs show less than 5% difference (see Figure 2.10), indicating that the 6-layer slab is thick enough to mimic GSFE in bulk aragonite. By comparing the NEB computed USFE structures with the perfect slab and the stacking fault structure, it is noteworthy that the USFE slab thickness increased by 0.22Å along z-direction, corresponding to a 0.73% expansion. To understand the effect of mechanical constraints, two more GSFE curves were calculated for the 6- layer slab as rigid sliding (no relaxation) and constrained sliding (only allowing layer expansion) as shown in Figure 2.10. The rigid sliding structures are built by simply changing the shift displacement between 0 and 4.161 Å upon perfect slab and 51 stacking fault slab without relaxation. The single point energies for these structures were calculated. The constrained sliding slab is constructed in the same way as the rigid sliding slab. But they are allowed to relax only along z-direction. The rigid sliding slab is built to mimic the deformation in mineral aragonite which does not have the nano-scale subgrain structure and biopolymers that allows the expansion along z-direction. The constrained sliding slab is mimicking the situation of no out-of-plane rotation of CO3 2- groups. Figure 2.10 The computed generalized stacking fault energy (GSFE) landscape and the fracture energy of aragonite. The USFE is higher than the fracture energy during rigid sliding, but drops below the fracture energy when layer expansion is allowed. Moreover, the USFE/SFE ratio is the highest in a broad range of metallic and oxide materials (as shown in Table 2.1). 52 Comparing the GSFE curves with and without constraints, the USFE dropped from 1400 to 900 mJ/m2 after allowing expansion along z-direction, and to 776.4 mJ/m2 after full atomic relaxation. In the USFE structure during rigid sliding, the CO3 2- is directly on top of each other and the closest O-O distance is only 2.04 Å, causing a large repulsion. If the layer expansion was allowed, this repulsion led to bending of the planar CO3 2- and increasing of the O-O distance to 2.66 Å in the NEB computed USFE structure. This causes the large energy drop of the USFE with layer expansion. More importantly, the computed fracture energy for aragonite is 980 mJ/m2. Therefore, if expansion along the [110] direction is not allowed (such as in mineral aragonite), a fracture will occur before sliding of (110) planes can be activated. However, the platelets in nacre are not a single crystal aragonite. Instead, it consists of nanocrystals with soft and viscoelastic biopolymer in between. As such, expansion along the [110] direction becomes possible in biomineral aragonite, and the USFE will be lower than the fracture energy. This renders activation of sliding in the form of stacking faults or partial dislocations prior to fracture. This mechanism may be the key difference between brittle mineral aragonite and toughened bio-mineral aragonite. Another bio- mineral, the bivalve Placuna placenta, also exhibits deformation twinning, while the mineral calcite is brittle[200]. For comparison, we also computed its SFE and TFE. The constrained USFE was obtained by the constrained method. 2.4.3. Compare the twinning tendency of CaCO3 with other materials The twinning tendency is a criterion to make comparisons between full dislocations (T>1) and twinning (T<1) based on GSFE curve as defined in Equation 2.8. It has been widely used to analyze the competition between full dislocations, stacking faults and twinning [201–205] in metals. It has been validated for many FCC metals, [206–208] 53 i=jk1+2∗l1(cid:28) (cid:1)‘mn(cid:1)o‘mnpq∗(cid:1)o‘mn (cid:1)ormn Equation 2.8 Tadmor and Bernstein[206] have calculated the twinning tendency for eight fcc metals. The result showed the metals have a twinning tendency order of Pt(cid:130)(cid:23)(cid:24)(cid:28)(cid:129)(cid:16)>(cid:130)(] where r is the interatomic distance and V is the potential energy. 12-6 Lennard-Jones potential is a pair-wise (or a two-body) potential that only considers the interaction between two atoms. These Equation 3.1 analytical models typically contain parameters that can be tuned to fit different systems. There are two parameters in 12-6 Lennard-Jones potential, ϵ, and σ. By calculating the derivative of energy with respect to atom position, the force between each atom pairs can be calculated. The total force on each atom is the summation of all pairwise forces. To increase the accuracy of simulation results, many forcefields have been proposed. M.P. Tosi et al. proposed Born-Mayer-Huggins (BMH) potential which considers only the two-body interactions in ionic systems [244]. This potential achieved relative accurate crystal lattice on various ionic materials with the error less than 0.05 Å and served as the basics for later developing of ionic system forcefield. BMH can also predict the silica structural reconstruction phenomena around the vacancy by the proper fitting of the forcefield parameters [245]. However, it is impossible for the two-body forcefield to have the correct description of bond angles, so three-body terms are added to ionic potentials such as the Vessal potential [246] which consists the Buckingham potential term and an additional angle term. 66 Feuston et al. [247] also proposed the Feuston-Garofalini (FG) potential which is a combination of BMH and Stillinger-Weber (SW) potential [248]. These three-body potentials can have correct predictions on the bond angle, radial distribution function, and structure factor comparing with diffraction and scattering results. In this study, we are implementing ReaxFF reactive forcefield which has the capability of modeling chemical reaction along with other materials properties. In ReaxFF, the total energy contains the following terms [131,249]: (cid:5)(cid:133)(cid:17)(cid:133)(cid:134)(cid:135)(cid:136)=(cid:5)g(cid:137)(cid:138)(cid:139)+(cid:5)(cid:137)(cid:140)(cid:135)(cid:141)+(cid:5)d(cid:138)(cid:139)(cid:135)(cid:141)+(cid:5)e(cid:142)+(cid:5)(cid:140)fe+(cid:5)(cid:134)(cid:137)(cid:141)+(cid:5)(cid:140)(cid:139)(cid:143)ffe(cid:133)+(cid:5)F(cid:137)de(cid:137)(cid:136)g Equation 3.2 where each term represents bonding energy, over-coordination penalty, under-coordination penalty, lone-pair energy, valence energy, torsion energy, van der Walls energy and Coulomb energy. In each interaction, the calculation is performed in the flow illustrated in Figure 3.2. All interaction terms can be classified into two groups: (covalent) bonded interactions and non-bonded interactions. For the bonded interactions, the first step is to determine the bond order from atomic positions. As shown in Figure 3.2, the bond order is determined as a function of interatomic distance, and the dependence on distance is continuous and smooth. To model the transition state, the bond order allows long term interactions. But this will lead to unrealistic interaction with the second neighbor of one atom. So the over-coordination correction term is introduced. The angle and torsion terms are three-body and four-body interactions based on the bonded pairs. So when a bond is broken, the bond order decreases, and the strength of a bond can decay continuously and smoothly. 67 Figure 3.2 Schematic of the ReaxFF iteration. The non-bonded interactions are on the left and the covalent/bonded interactions are on the right [131] The non-bonded interactions including Coulomb interaction and van der Walls interaction. These two energies are calculated across all atom pairs regardless of their bonding situation. They are also shielded at the short distance by introducing a shielding term to avoid unphysical large values due to the small distance on the denominator. The atomic charges are dynamically allocated by the electronegativity equalization method (EEM). More details of the EEM method will be discussed in Chapter 5. 68 Figure 3.3 The bond order of a carbon-carbon interaction with respect to distance. [131] Based on the analytical formulation of the forcefield, the potential energy of the system and the forces on each atom can be determined. The dynamic motion of the atomistic system will be tracked over time. For the given atomic system, the temperature is defined as the kinetic energy of a system by Equation 3.3: T= 23Nk(cid:148)(cid:5)U Equation 3.3 69 where N is the number of particles, ES is the total kinetic energy and k(cid:148) is the Boltzmann constant. At the beginning of dynamics, all atoms are assigned an initial velocity randomly based on the desired simulation temperature. The motion of the atoms is governed by the Newton’s Law. To actually solve the trajectory of each atom as different time steps, the problem becomes an initial value problem. The numerical solution is computed by Verlet algorithm. For the microcanonical ensemble, the total energy, total number of atoms, and volume of the system are constant. Thus it is also called NVE ensemble. Usually, the temperature will deviate from the initial temperature due to the transfer of energy between kinetic and potential energy. In order to achieve a relative constant temperature during the simulation, NVT (constant atom number, volume, and temperature, also refer as canonical ensemble) is applied. Instead of simply scale all velocities together according to the desired temperature and result into unphysical results, a thermostat is used to achieve the trajectory that is consistent with statistical mechanics of the NVT ensemble during adjusting the velocities. Currently, Nosé–Hoover thermostat is widely implemented in the simulation to control the temperature[92,250]. In Nosé–Hoover thermostat, fictitious thermal bath is added to the existing system, and a new Hamiltonian is defined for the extended system.. The constant pressure can be achieved using NPT ensemble (constant atom number, pressure, and temperature, also refer as isothermal-isobaric ensemble) by applying a barostat using the similar idea as thermostat. 3.3.2. Simulation of the forming process of oxide bifilms In this chapter, all MD simulations were carried out using LAMMPS package[251,252]. The parameters used in this study were presented and used by Sen et al. [134] to investigate the oxidation of crystalline aluminum nanowire. 70 Figure 3.4 The procedure to build aluminum oxide bifilm structures with different formation history To compare the change in oxide bifilm mechanical properties during the aging process, a protocol, shown in Figure 3.4, was developed for the aluminum oxide bifilm structures with different formation history. The corresponding structure details are listed in Table 1. As illustrated in Figure 3.1, the formation history of the bifilm is complicated in terms of both time and temperature change. However, the typical MD timescale (~ns) is too short to observe the crystallization and the aging 71 process of aluminum oxide occurring in hours. Therefore, two representative structures and conditions, young oxides and old oxides, were built on aluminum metal. Table 3.1 Simulation details for monolayer and oxide bifilm structures shown in Figure 3.4 Number of Number of Number of Structure name Figure 3.4(c) Monolayer young oxides Figure 3.4(d) Monolayer old oxides Figure 3.4(e) -OH Monolayer old oxides terminated Figure 3.4(f) Young oxides Figure 3.4(f) Old oxides Figure 3.4(f) -OH terminated old oxide Al 6708 O 1539 7885 1879 7885 13416 15770 15770 1879 3077 3758 3758 H - - Cell size/Å 80×40×140 80×40×140 200 80×40×140 - - 78.49×39.93×140 81.18×41.64×140 400 80.86×41.49×140 As shown in Figure 3.4(a), for young oxide, its formation was directly simulated using MD simulation by allowing an 8 nm thick Al (1 0 0) slab reacting with the O2 environment. To accelerate the oxidation process, the pressure of O2 was set to 200 atm and the temperature was set to the melting point of Al. After 200 ps NVT (constant number of particles, volume, and temperature) MD simulations, a ~1 nm thick aluminum oxide passivating layer was formed on the Al surfaces on both sides. Then, a 4 nm thick slab (3 nm thick of aluminum covered by 1nm aluminum oxide) was taken from the oxidized Al-slab as the structure of monolayer young oxides (see Figure 3.4(c)). 72 As shown in Figure 3.4(b), the old oxide structure was built by attaching 1 nm thick #(cid:28)hw(cid:24)x\ to a 3nm thick Al (1 0 0) slab, then relaxing the structure with 50 ps NVT MD simulations at the melting point of Al. However, the sharp interface between the Al and the oxide indicates limited inter-diffusion at the melting point of Al (See section 3.2 for more details). To have a fully smeared Al/oxide interface that is formed due to oxidation, the Al/oxide interface was heated to the melting point of aluminum oxide, as indicated in Figure 3.4(d). The –OH terminated old oxide surface was simulated by terminating ~30% surface oxygen atoms with hydrogen atoms. Therefore, 200 hydrogen atoms were added to the old oxides monolayer slab surface. The –OH terminated old oxides structure was minimized and equilibrated using NPT (constant number of particles, pressure, and temperature ensemble to release any stress in the x-y plane) ensemble at the melting temperature of Al. The relaxed structure is shown in Figure 3.4(e). Finally, three monolayer oxide/Al slabs were formed via a series of MD simulations and the Al/oxide interface in all three cases was considered to be the “wetted” side of the oxide[221]. The folding process that forms the oxide bifilm was simulated by duplicating the monolayer, flipping it and placing it on top of the original monolayer slab. Forces that correspond to 1 GPa pressure was applied to the aluminum atoms at the top and bottom of the bifilm slab for 50 ps using NVT ensemble followed by 10 ps NPT (to relax the stress on the x-y plane). As discussed earlier, more bonding may form across the bifilm interface with increasing aging time and temperature. Therefore, to compare the two extreme folding conditions, young oxides were folded at the room temperature (and compared with folding at high temperature) and old oxides were folded at the melting temperature of Al then quenched to the room temperature followed by a 10 ps NVT MD relaxation (see section 3.2 for more discussion). The –OH terminated old oxide was folded at room temperature because the H2 will only escape from aluminum matrix after 73 solidification at low temperature. The interface between the folded aluminum oxide slabs is considered to be the “dry” or “inner” interface[221]. 3.3.3. Simulation of the fracture process of the aluminum oxide bifilms Uniaxial tensile simulations at room temperature were performed along the direction that is perpendicular to the slab surface. The length of the simulation cell was increased by 0.5% of strain along the tensile direction every 1ps, followed by structure relaxation at room temperature with NVT dynamics. To maintain the applied strain in a cell with vacuum, the atoms in the top and the bottom layers in the bifilm structure were fixed in the z-direction during structure relaxation. This simulates an applied constant tensile strain rate of 5×109/s. 3.4. Validation of ReaxFF To validate the accuracy of this ReaxFF reactive force field [134], the modulus, melting temperature and density of different oxide structures were calculated and compared with experimental values. The results are described in this Section. 3.4.1. Scaled temperature The melting temperature of Al was first calibrated in this study. To avoid the overheating and overcooling effect in MD simulation due to the lack of defect as nuclei, we implemented a two- phase method to calculate the melting point[253], which is defined as the temperature for liquid and solid to co-exist. In the two-phase method, the crystal Al and liquid Al was built separately at 700K and 800K. Then the liquid structure was quenched to 400K and stacked with crystal Al to form a solid-liquid interface. The different initial temperature was applied to the two-phase structure with NVE (constant number of Al, constant volume and constant total energy) ensemble. 74 The interface will move towards the solid side or the liquid side depending on the different initial temperature. When the final temperature reached equilibrium with a clear solid-liquid interface, this final temperature is considered to be the melting point of Al. The details of two-phase simulation and a case with low initial temperature is shown in Figure 3.5. It clearly showed the liquid phase solidified and the crystalline phase grew and the solid/liquid interface move toward the liquid phase. Several stacking faults (red atoms) were formed during crystallization. The melting point of Al calculated using the ReaxFF is 670K, which is lower than the 933K experimental melting temperature. This is partly because the melting point data is not explicitly included in the ReaxFF training set. To compensate for this effect, we used scaled temperature in all simulations in this study. A similar approach has been taken in previous MD simulations[59,134]. Considering the room temperature, 300K, which is about 32% of actual Al melting temperature, we used 32% of calculated Al melting temperature as the scaled room temperature, meaning RT=210 K=32% Tm(Al) ≈ 32% × 670 K. 75 Figure 3.5 Build structure for two-phase simulation. The initial temperature of 200K is applied to show the moving of interface toward crystal side 3.4.2. Density and Elastic Modulus of the oxides The density is validated across different aluminum oxide phases. The density is calculated by applying NPT dynamics to bulk Al(cid:24)O\ structure until the system reaches equilibrium. At the scaled room temperature, the calculated densities of amorphous, γ(cid:28) and α(cid:28)Al(cid:24)O\ are 3.47, 3.46 and 3.74 g/cm3, respectively. The experimental values are 3-3.5, 3.66 and 4 g/cm3, 76 respectively[254,255]. The relative errors are less than 7%. At the scaled melting point of 3.72 g/cm3 , which are reasonably smaller than the density at room temperature due to volume aluminum, the corresponding density values of amorphous, γ(cid:28) and α(cid:28)Al(cid:24)O\ are 3.45, 3.43 and expansion. The corresponding linear thermal expansion coefficients are 4.2×10-6/K, 6.3×10-6/K, and 3.9×10-6/K, which are consistent with the experimental value of 7.5×10-6/K[256]. Tensile simulations were performed for amorphous-Al2O3 at scaled room temperature (0.3 Tm) and Al melting temperatures (Tm) to obtain the Young’s modulus. According to the predicted Tm, room temperature (0.3 Tm) is 210K. Amorphous structure with 5500 atoms of Al2O3 is generated by heating κ-Al2O3 crystal structure to 4000K then rapidly quenching to the 0.3 Tm and Tm. For comparison, a series of structures with a density ranging from 2.0 g/cm2 to 4.0 g/cm2 at different O/Al ratios ranging from 1.1 to 1.5 is built by resizing the cell length and randomly deleting O atoms. Before deformation, each structure is heated to 4000K for 10 ps to form the liquid and quenched down to 0.3 Tm and Tm for 10 ps to reach equilibrium with NVT dynamics. The tensile simulation is performed along z-direction while the cell lengths in the transverse direction are fixed. The relation between modulus and Poisson's ratio for an isotropic material is described in Equation 3.4: (cid:151)(cid:152)(cid:23)(cid:23)(cid:152)(cid:24)(cid:24)(cid:152)\\(cid:153)=1(cid:5) (cid:151)1(cid:28)(cid:154)(cid:28)(cid:154) (cid:28)(cid:154)1(cid:28)(cid:154) (cid:28)(cid:154)(cid:28)(cid:154)1(cid:153)(cid:151)(cid:16)(cid:23)(cid:23)(cid:16)(cid:24)(cid:24)(cid:16)\\(cid:153) Equation 3.4 where σ, ε, ν and E are stress, strain, Poisson’s ratio, and Young’s modulus, respectively. Since we have applied (0, 0, (cid:152)\\) strain and computed the stress ((cid:16)(cid:23)(cid:23),(cid:16)(cid:24)(cid:24),(cid:16)\\) with NVT dynamics, the Young’s modulus can be obtained by fitting the stress-strain data and the result is shown in Figure 77 3.6. In the experiment, the Young’s modulus for amorphous Al2O3 with 2.8 g/cm2 density is ranging from 95-110 GPa [257,258]. Using NPT dynamic, the equilibrium density of Al2O3 at room temperature is about 80 GPa. Comparing with the experimental data, the ReaxFF predicted about 10% lower result, which is within a reasonable deviation for modulus calculation. Figure 3.6 The ReaxFF predicted room temperature Young's modulus of Al oxide with different oxygen content at different density The elastic bulk modulus of the three aluminum oxides was calculated by fitting the energy of the bulk structure under hydrostatic deformation in the range of (cid:155)5% strain. These fully periodic structures, without any defects, tend to deform mainly elastically within this strain range. The calculated bulk modulus of amorphous, gamma and alpha aluminum oxide are 80, 143 and 191 GPa, as plotted in Figure 3.7. The open and solid dots are moduli data measured by experiments or calculated using DFT method from the literature [233,259–270]. Although the absolute modulus values are scattered, the order of the modulus in amorphous, γ− and then α− (from low to high) is consistent. The absolute value of the bulk modulus calculated in this work is slightly lower than 78 the literature data but the trend reminds the same. Overall, the ReaxFF can provide reasonably accurate quantitative results on the bulk modulus. Figure 3.7 The predicted room temperature bulk modulus of bulk aluminum oxides at different phases in comparison with experimental values in the literature [233,259–270] 79 3.5. Results and Discussion 3.5.1. Interface analysis in oxide bifilms To characterize the “wetted” Al/oxides interface and the “dry” oxide/oxide interfaces, the atomic charge distribution and the bond density at these interfaces were analyzed for the fully-developed aluminum oxide bifilm structures. Figure 3.8 shows the atomistic structures of the “wetted” Al/oxide interface in the mono-slab and the projected charge distributions normal to the interface. The z position, or the surface normal direction, is labeled on the vertical axis of Figure 3.8. The right panel shows the relaxed interface structure with red, gray and black points representing oxygen, aluminum and hydrogen atoms. Their “atomic charges” are projected on the left panel with the same color scheme. The oxygen charges are about -1.0 in the oxide layer and the charges on Al atoms change from neutral in the metallic phase to as high as +1.5 in the oxide phase. The Al/oxide interface region has a finite thickness, which is bracketed by the dashed lines highlighting the positions of the top of the neutral Al atoms and the lowest position of oxygen in aluminum. For clarification, the structures in Figure 3.8(a), (c) and (d) are room temperature ones and the structure in Figure 3.8(b) is at the melting temperature of Al. Figure 3.8(a) shows a ~1.1nm thick amorphous young oxide layer on the top of an aluminum slab, under which there is an Al/oxide interface region with a thickness of ~0.8 nm. This interface region was formed during direct oxidation of the aluminum metal in an oxygen atmosphere at the melting point of Al. The oxidation process was accompanied by the inter-diffusion of oxygen and aluminum atoms, which generated a smeared and well-bonded interface. 80 Figure 3.8 The charge on each atom of (a) the young oxides, (b) old oxides without extra heating, (c) old oxides and (d) –OH terminated old oxides with respect of their position along the surface normal direction. The corresponding atomic structures are shown 81 Figure 3.8(b) shows the mono-slab structure after relaxing the as-built amorphous Al2O3/Al interface for 50 ps at the melting temperature of aluminum. Comparing to the monoslab with young oxide, the interface thickness region in Figure 3.8(b) is only 0.2 nm. This sharp interface is due to insufficient interdiffusion within the limitation of total simulation time in MD simulations even at the melting temperature of aluminum. Considering the actual interface of the old oxides grown on aluminum should have at least similar smeared thickness with the young oxides, if not thicker due to its longer formation time, a local heating process was added to enhance the interdiffusion at the Al/oxides interface (step (d) in Figure 3.4). More specifically, a 7 Å thick interface region was heated to the melting temperature of aluminum oxide (Al2O3) for 50 ps while the bulk Al part is still at room temperature. This resulted in the old-oxide/Al interface structure shown in Figure 3.8 (c), in which the interface region thickness increased to 0.6 nm. This thickness is comparable with that in young oxides formed due to oxidation. Figure 3.8(d) shows the –OH terminated oxide/Al interface at room temperature. The H atoms are ~ +0.4 charged and located within the top layer of the oxide. It has a width in the charge distribution plot because the top surface is not smooth. The top position of the oxide surface, in fact, fluctuates within ~5Å. Since it was built upon the relaxed interface structure in Figure 3.8 (c), the interface region thickness remains 0.6 nm. So, all three structures have well “wetted” and smeared Al/oxides interfaces. 82 Figure 3.9 (a) The Al-O bond density of (b) young oxides, (c) old oxides and (d) -OH terminated old oxides and their bifilm structures The mono-slab structures shown in Figure 3.8 (a), (c), and (d) representing the young oxides, old oxides, and –OH terminated old oxides, are then folded to form the oxide bifilms. Figure 3.9 shows the fully relaxed bifilm structures and the Al-O bond density distribution normal to the interface for these three types of oxide bifilms. The bond density is the number of Al-O bonds (with a cutoff 83 distance of 2 Å) in the cross-section parallel to interface plane. In Figure 3.9 (a), the blue, red and yellow curves indicate the Al-O bond density of young oxides, old oxides and –OH terminated old oxides at the dashed box regions, respectively. The zero point on the vertical axis corresponds to the middle plane of the oxide/oxide interface. As a comparison, the bond density in the amorphous aluminum oxide is plotted in dash line. In the Al layer, only the Al-Al bond exists. So, the bonding density of Al-O decreases to zero at the two ends. At the Al/oxide interface, Al-O bond density is increased gradually. The gradual transition is also consistent with the smeared Al/oxide interface. The dashed line indicates average bond density for bulk amorphous aluminum oxide. It is slightly lower than the bond density at the center of the oxide layer since the oxide slab is compressed in the folding process. It also indicates that the oxide slab thickness is sufficient, as its structure can represent the bulk structure. The “dry” oxide/oxide interface has much lower but non-zero Al-O bonds density. This means that the oxides are partially “healed” by forming some Al-O cross the “dry” interface for all three structures. The bond density of young oxide is lower than that of old oxides. This is because young oxides are folded at a lower temperature. We also simulated folding at the melting temperature of Al. By folding young oxides at high temperature, the bond density increased to 0.06/A2, as the same as the old oxides. Also, the oxygen content difference on the “dry” interface contributes to the lower bond density. During the oxidation process, the oxygen content in the surface oxides increased gradually. For young oxides, the actual O/Al ratio is less than 1.5, which was also observed in other MD simulation studies[134] and experiments.[271] But the old oxide retains the exact 2:3 Al/O ratio as it is the α-Al2O3 and more Al-O bonds can form at the interface. The –OH terminated old oxide interface has a much lower bond density comparing with the old oxide 84 interface due to the hydrogen atoms on the surface preventing the Al-O bond formation across the “dry” interface. In all three oxide bifilm structures, the Al-O bond density at the “dry” interface is about 50% to 75% of the bond density in the bulk amorphous Al2O3. This means that the two oxide layers will not fully “heal” although a pressure of 1 GPa is applied during the bifilm formation process. Even though there are a partial “healing” process and some Al-O bond forms across the oxide/oxide interface, the oxide bifilm indeed has a two-layer structure with significant low bonding density gap in the middle. This is consistent with the “non-wetted” assumption of the inner interface of the bifilms.[221,272,273] 3.5.2. Deformation of the aluminum oxide bifilms under tension The oxide bifilm structures, including the two kinds of interfaces in the whole Al-oxide-oxide-Al slab, were subjected to uniaxial tension at room temperature. Figure 3.10 (a) shows the fracture feature of young, old, and -OH terminated old oxide bifilms at 20% strain, where a fully detached crack surface can be observed. Note these oxides are much more ductile than their bulk form or even considered to be “super plastic”[59,134] or “liquid-like”[60], therefore a large strain is required. It is seen that both young and old oxide bifilms fracture at the “wetted” Al/oxides interface, and some aluminum metal transfers to the oxide after fracture. Unlike young and old oxides, the –OH terminated old oxides fracture at the “dry” oxide/oxide interface. 85 Figure 3.10 Deformation analysis of bifilms. (a) The fracture feature of young, old and –OH terminated old oxide bifilms. (b) The stress-strain curves of the three bifilms. The stress-strain curves of the three oxide bifilms are plotted in Figure 3.10 (b). The blue, red and yellow curves represent young, old and –OH terminated old oxides, respectively. During the deformation process, the stress for the old oxide bifilm increases to the ultimate tensile stress and then drops immediately to zero stress at ~6.5% strain, exhibiting a brittle-fracture behavior. The 86 following stress oscillation is due to the stretching of several aluminum nano-ligaments developed during the fracture process, as shown in Figure 3.10(a). For the young oxide bifilm, a small stress plateau following the ultimate tensile stress indicates a more ductile tensile behavior comparing with the old oxide bifilm. However, the simulated ultimate tensile stress of young oxides is 1.1 GPa, which is 0.45 GPa lower than that of the old oxides (1.55 GPa). For the –OH terminated oxide bifilm, the ultimate tensile stress is 1.2 GPa, similar to that of the young oxide bifilm, but the stress drops gradually, due to the plasticity exhibited in the nano-scale oxide films. Similar to young and old oxides, several nano-ligaments between the fracture surfaces are still connected after fracture. For all three oxide bifilms, major cracks are formed across the fractured interface around 6% strain. We notice that although the Al-O bond density at the oxide-oxide “dry” interface in the bifilm is lower than the rest of the oxides, the oxide/oxide interface is still stronger than the Al/oxides interface. Increasing the folding temperature of the young oxide will not change the location of the fracture, as it will further increase the bond density and enhancing the bonding at the “dry” interface. This is different from the general belief that the “wetted” Al/Al2O3 interface is well- bonded and stronger than the “dry” side. In fact, the strength of metallic Al-Al bonds is intrinsically weaker than the ionic Al-O bonds, supported by a density functional theory (DFT) study by Wang et al.[138] They compared the work of separation at different positions in the Al/oxide slab and found that the Al/oxide interface had the lowest work of separation of 0.79 J/m2. But the work of separation only characterizes the energy change in equilibrium state since the DFT calculations did not take the temperature effect, nor plastic deformation into consideration. In the MD simulation, the fracture energy can be calculated by integrating the area under the stress- strain curve from zero strain to the fracture strain. In this way, the temperature effect and energy 87 change due to plastic deformation can be included. In our study, the calculated fracture energy for young and old oxide slabs is 0.43 and 0.53 and J/m2, respectively. The fracture energy, especially for the one fractures at the Al/old-oxide interface, is in good agreement with the DFT results. The fracture energy of the –OH terminated old oxide is only 0.30 mJ/m2, which explains why fracture occurs at the “dry” oxide interface. Also, all three structures show higher ductility than typical experimental observation of bulk aluminum oxide. This is because, at the nanoscale, the surface diffusivity is greatly enhanced. During the deformation process, the Al-O bond can be stretched to break and then healed by reconnecting with neighboring atoms.[134] This is similar to the findings reported for the compression test on aluminum oxide nanoparticle[241,242] and nano-indentation on TiO2.[274] Comparing our simulation results with experimental observations also reveals some interesting conclusions. Wang reported that after deformation of the casting samples in the experiment, young oxides always appear on one side of the fracture surface, which is consistent with our simulation observation. But old oxides always appear on both sides of fracture surfaces which means old oxides fracture at the oxides/oxides interface.[147] Also, SEM images acquired by Dr. Jianfei Sun from Harbin Institute of Technology on fracture morphology of bifilms at different stage agrees well with our MD simulations. Our calculation shows that the phase transition from amorphous oxide to α-Al2O3 and the folding temperature from the room temperature to the melting point of aluminum will not change the location of the fracture. Therefore this difference is more likely due to the impurities or the hydrogen trapped at the old oxide bifilm interface during the aging process. It has been suggested in the literature that the –OH termination on the contaminated surface repels each other and at the same time occupies possible reactive sites for Al-O bond formation during the folding 88 process.[236–238,275,276] This is proved by the Al-O bond density plot shown in Figure 3.9. The H atoms at the bifilm interface also prevent dangling Al-O bonds (after bond breaking) to be reconnected with other Al and O atoms,[134] thus reducing the fracture energy. Overall, the fracture energy of the –OH terminated old oxides is much lower than that of the non-contaminated old oxides. Therefore, the simulation results further conclude that the contribution of –OH group surface contamination to the initiation of cracks at the “dry” interface is more significant than the reduced Al-O bond density and the brittle nature of the oxide. 3.5.3. Scale bridging to microscale for cohesive zone model Simulations of the oxide distribution during the casting process and its fracture in the solidified cast aluminum parts are both important. Pita and Felicelli have proposed an immersed element- free Galerkin (IEFG) method to model the deformation of the oxide film in liquid aluminum[150]. They assumed the total volume and density of oxide bifilms remain unchanged due to the lack of information on oxide properties. In most engineering applications, finite element method (FEM) is used for solid structure analysis and a cohesive zone model (CZM) can be adapted for deformation and fracture calculation.[277] Therefore, it will be extremely useful to pass the MD predicted aluminum oxide bifilm fracture properties to the CZM and other continuum models for engineering design of aluminum castings. It is noted that the predicted fracture strengths of all three aluminum oxide structures are above 1 GPa which appears to be much higher than the measured values for a typical metal/oxide interface (hundreds of MPa). This is a well-known size effect due to two reasons: first, the lack of defects such as stacking faults and dislocations to initiate the yield process[278]; second, the limited size of simulation cell at atomistic level forbidding the gradual propagation of cracks. In the FEM model, the smallest step for crack propagation is the mesh size. In MD simulation, due to the cell size and the periodic boundary condition, the crack 89 tends to open up the whole interface at the same time. Xia et al. compared the MD predicted interface property and the CZM parameters needed to predict the experimental measurements for the Al/Si interface in cast Al-Si alloys.[279] They showed that the fracture energy is conserved in both MD and CZM simulations. They proposed the relationship between fracture strength in the MD simulation, (cid:16)∗, and the fracture strength in FEM, (cid:16)(cid:19), can be scaled by Equation 3.5: lnl(cid:158)(cid:159)p=(cid:28)12ln=(cid:16)(cid:19)(cid:16)∗G Equation 3.5 where L is the mesh size in the FEM model and a is the intrinsic cohesive zone of length. Figure 3.11 The fracture strengths of the three bifilms under different mesh size in finite element methods. 90 Assuming the oxide bifilm will be treated as crack nucleation in FEM with a CZM [279], we scaled the fracture strength with the mesh-size used in FEM, as shown in Figure 3.11. More theory derivations for this size bridging method can be found in Ref. [279]. The results showed that under 100nm mesh size, the fracture strength of young, old and –OH terminated old oxides are 50.7, 71.1 and 63.7 MPa, respectively. These values can be used as basic oxide properties to predict the fracture behavior of aluminum castings. Due to the difficulty of conducting experiments to probe the mechanical properties of the oxide bifilms [280], there is no direct verification of our prediction yet. This is indeed the opportunity and motivation to apply accurate atomistic modeling to predict properties that are hard to measure. The MD simulation results are consistent with several previous experimental and computational studies indirectly. For example, Nyahumwa [281] reported that in unfiltered casting samples (consist of both young and old oxides), fractures occur at the young oxides. This indicated that the fracture strength of the young oxides should be lower than that of the old oxides, which agrees with our MD results. In the IEFG model [150], it was assumed that the oxide bifilm was stretchable, which is consistent with the ductile behavior observed in MD simulations. Based on equation 1, it can be derived that when mesh size is at the magnitude of a millimeter, the corresponding fracture energy is on the order of 5×10(cid:22)(cid:160)MPa, which is close to zero. This is the assumption used for analytical fatigue life models, which assume that fatigue crack initiation life is zero if the crack is initiated from defects like porosity and oxides[239,282]. These agreements collectively supported the simulation results. Therefore these MD results can be used to quantify oxide bifilm deformation and fracture properties at different aging stage and different chemical environments. These predictions will then serve as inputs to a multiscale modeling framework to simulate the oxides distributions in the casting process and the fracture of the casted components during operation. 91 3.6. Conclusion Based on the formation history of oxide bifilms, a set of reactive MD simulation procedures were carefully designed to predict the mechanical properties of the oxide bifilms formed in aluminum casting. During oxide bifilm formation, the Al-O bond density at the “dry” interface is about 50% to 75% of the bond density in the bulk amorphous Al2O3, indicating an incomplete “healing” process occurred in the oxide. During the oxide aging process, the oxide transforms from amorphous to α-Al2O3 and the fracture energy also increases from 0.43 J/m2 to 0.53 J/m2. Despite the “healing” process, both young and old oxide bifilms fracture at the “wetted” Al/oxides interface with some aluminum metal transferred to the oxide after fracture. In comparison, with 30% coverage of hydroxyl group surface contamination and similar Al-O bond density across the oxide/oxide interface, the –OH terminated (old) oxide bifilm tends to fracture at oxide/oxide interface and the corresponding fracture energy drops to 0.30 J/m2. The MD predicted fracture strengths are 1.1 GPa, 1.55 GPa, and 1.2 GPa for young, old, and –OH terminated old oxide bifilms, respectively. The hydrogen is likely to come from the H2 bubble trapped in the aluminum oxide bifilm interface. Therefore, the simulation results indicate that the contribution of –OH group surface contamination to the initiation of cracks at the “dry” interface is more significant than the phase change during aging, the interface healing process during bifilm formation, and the brittle nature of the oxide. Therefore, removing H2 trapping, via degassing or adding compounds that react with hydrogen more favorably than the oxide, will reduce a possible source for –OH termination, thus allow the oxides heal and strengthen at the bifilm interface. For large scale cast aluminum part design, the oxide bifilm initiated crack opening can be described by a cohesive zone model (CZM) in finite 92 element method (FEM). Therefore, a simple size-bridging relationship was used to scale the MD predicted fracture strengths to CZM parameters. 93 CHAPTER 4. Oxidation modulated Lithium nano-structure growth morphology change during extrusion (Part of this chapter is adopted from Ref. [151]) 4.1. Summary Li metal is the preferred anode material for all-solid-state Li batteries. However, a stable plating and stripping of Li metal at the anode-solid electrolyte interface remains a significant challenge, particularly at practically feasible current densities. This problem is usually related to high and inhomogeneous Li-electrode-electrolyte interfacial impedance and to formation and growth of high aspect ratio dendritic Li deposits at the electrode-electrolyte interface which eventually shunt the battery. To assist the understanding the details of Li metal plating, operando electron microscopy and Auger spectroscopy were used to probe nucleation, growth, and stripping of Li metal during cycling of a model solid-state Li battery as a function of current density and oxygen pressure. Molecular dynamics simulations were performed to illustrate the impact of the oxygen partial pressure on the Li-growth morphology. We find a linear correlation between the nucleation density of Li clusters and the charging rate (growth rate) in an ultra-high vacuum, which agrees with a classical nucleation and growth model. However, the trace amount of oxidizing gas (≈ 10-6 Pa of O2), promotes the Li growth in the form of nanowires due to fine balance between the ion current density and growth rate of a thin lithium oxide shell on the surface of the metallic Li. Interestingly, increasing the partial pressure of O2 further (to 10-5 Pa) resumes Li plating in the form of 3D particles. Our results demonstrate the importance of trace amounts of preexisting or ambient oxidizing species on lithiation processes in solid-state batteries. 94 4.2. Introduction and Experimental Observation 4.2.1. The Li metal and Li2O in battery applications Li metal is an attractive anode material for all-solid-state Li batteries (SSLBs) due to its high theoretical capacity (3.86 mAhg-1) and low potential vs. standard hydrogen electrode (-3.04 V)[283,284]. In order to match the practical energy capacity, a relatively thick Li metal film of ≈ 20 µm is needed for a typical Li-ion SSLB[285]. The current major experimental effort is dedicated to achieving homogeneous plating/stripping of such a thick Li layer to avoid destructive mechanical strains in the anode and to minimize the impedance of electrode/SSE interfaces.[286– 290] Though significantly mitigated compared to liquid electrolytes, the shunting of the battery by metallic Li filaments growing through grain boundaries of crystalline SSEs remains to be an issue.[291] During the battery charging and discharging process, one of the most severe safety concern is the dendrite growth problem at the Li/SSE. [292–294] After serval cycles, the Li dendrite formed on this interface will penetrate through the electrolyte and cause short circuit. The aforementioned SSLBs challenges become particularly important for current densities exceeding ≈ 1 mA/cm2, which is still well below the >10 mA/cm2 requirement for portable electronics and electric vehicle applications.[291] Another challenge is the capacity loss due to the side reactions that form the solid electrolyte interphase (SEI). This thin layer should be electronically insulating and Li-ion conductive. Due to the complex morphology change of Li anode during cycles, achieving a stable plating and stripping of Li metal deposits in contact with a solid-state electrolyte (SSE) in such a cell yet remains a significant technological challenge, even when the Li/SSE interface is thermodynamically or kinetically stabilized [285,295]. 95 The lithium oxide is a major component found in all the solid-electrolyte interphase (SEI) for many energy storage devices such as for Li-ion battery [296], all-solid-state battery, and Li-air battery [297]. Due to the high chemical reactivity and oxygen affinity of Lithium metal, the oxidation of lithium is inevitable [57,58] and result in nanometer thick lithium oxide. Due to the volume change of the lithium electrode, the previously formed passivating SEI layer may still crack during the charging and discharging process, causing degradations of battery performance [298–300]. Thus, it is important to first use Li2O as a repressive SEI layer and investigate its formation and impact on Li-morphology change dynamically. The coupled chemical-mechanical process will be insightful for battery design. Figure 4.1 (Left) A schematic of the experimental setup used to cycle the all-solid-state batteries in an operando SEM with the controlled pressure of O2. The grounded tungsten (W) tip is used to contact the top C-Al anode. A potential is applied to the bottom Pt-Ti current collector electrode to charge/discharge the electrochemical cell. (Right) Real-color optical image of C-Al anode in a representative device. 96 4.2.2. Experimental observations To further understand the evolution of Li anode morphology during the charging process, experimental techniques are developed to characterize Li metal plating and stripping on a model ultra-thin carbon anode deposited on amorphous, nitrogen doped lithium phosphate (LiPON) SSE /LiCoO2 cathode stack. A schematic of the thin film battery and the experimental set-up are illustrated in Figure 4.1. Arrays of SSLBs are fabricated on a Si (001) wafer covered with a 100 nm thermal SiO2. A 90 nm Pt /30 nm Ti current collector is deposited using an electron beam evaporation followed by a ≈345 nm thick LiCoO2 cathode sputtered using previously described conditions.[301] Following deposition, the LiCoO2 is annealed in the O2 atmosphere at 700 oC for 2 h to form the high-temperature phase. A ≈340 nm thick LiPON electrolyte layer is sputtered on top of the cathode layer. Next, an array of amorphous carbon (C) anodes with a thickness of ≈ 35 nm is evaporated through a stencil mask with a diameter of 0.51 mm. Then ≈ 210 nm thick Al pads have been evaporated through the same shadow mask but with an offset to leave most of the C anode surface exposed. The Al pads are used to make electrical contact with the C anodes. To elucidate the role of oxidizing species from the ambient environment, The Li morphology change at different O2 partial pressure in the range from ≈ 5 × 10-7 Pa to ≈ 5 × 10-5 Pa is studied. The observations showed that the morphology of metallic Li deposition is highly sensitive to the oxidizing ambient and drastically changes from mostly in-plane 3D particle growth mode under UHV conditions to out-of-plane whisker-like growth once the oxygen pressure in UHV chamber increases from ≈ 10-7 Pa to ≈ 10-6 Pa. The in-plane Li growth mode reestablishes once the oxygen pressure further rises to 10-5 Pa. 97 SEM images in Figure 4.2 depict the galvanotactic Li plating at 0.77 mA/cm2 current density under 5.7 × 10-7 Pa, 5.7 × 10-6 Pa, and 4.8 × 10-5 Pa oxygen pressure, respectively. Interestingly, Li deposit morphology changes from the mainly in-plane 3D growth of metallic Li islands with a negligible fraction of nanowires, to one dominated by nanowire growth as the O2 pressure increases one order of magnitude from ≈ 6 × 10-7 Pa to ≈ 6 × 10-6 Pa. In-plane Li 3D particles growth mode is again re-established as the O2 partial pressure is further increased to ≈ 5 × 10-5 Pa. Moreover, the nucleation density of Li particles increases three-fold across the oxygen pressure span from ≈ 6 × 10-7 Pa to ≈ 5× 10-5 Pa (Figure 4.1 a-c). Figure 4.2 Effect of oxygen pressure on Li plating morphology. SEM images of Li plated at (a) 5.7 × 10–7, (b) 5.7 × 10–6, and (c) 4.8 × 10–5 Pa residual oxygen pressure, respectively; 0.77 mA/cm2 current density. The SEM images are acquired at ≈ 0.04 mAh/cm2 capacity of deposition Li nanowire growth mechanism and critical thickness of the oxide layer It is postulated that the change in the Li growth mode to the formation of an encapsulating lithium oxide surface layer that controls the morphology of growing Li nanoparticle. But due to the extremely thin thickness of surface oxides and small diameter of Li nanowire, it is challenging to conduct in-situ experiment to observe and quantify the growth mechanism of Li nanowire during 98 the charging process. Therefore, to understand the competition between oxidation rate and current density, ReaxFF reactive forcefield based MD simulation is implemented in this study. The methods to control oxidation rate and Li growth rate is introduced in Chapter 4.3. Then the observation in MD simulation is presented in Chapter 4.4.1. Based on MD observation, the Li growth mechanism and reason for its morphology change is discussed in Chapter 4.4.2. Also, we predicted the morphology of Li nanowire change under different experimental conditions in Chapter 4.4.4. 4.3. Computational details To simulate the growth of Li nanowire under different O2 partial pressure and current density, ReaxFF reactive forcefield based molecular dynamics is implemented[302]. Different from classical forcefield, the ReaxFF reactive forcefield can predict the charge on each atom and model the chemical reaction between different molecules. 4.3.1. Oxidation of Li and Al slab In the experiment, Li metal will be rapidly oxidized when exposed to air or even under ultra-high vacuum condition. So, the surface is also oxidized in MD simulation for 10 ps at 10 atm under 150 K to mimic the thin oxide layer formed in the experiment. In our study, an 8nm×8nm×9nm Li slab was built first. Then, the slab is placed in a simulation box with the size of 8nm×8nm×20nm as shown in Figure 4.3. To simulate the experiment environment, a certain number of oxygen molecules were randomly distributed in the vacuum region to achieve a given pressure. For example, 335 O2 molecules in the volume of 8nm×8nm×10nm = 6.4×103 nm3 empty space in the simulation cell to achieve 10 atm at 150 K. The high pressure is used to achieve a reasonable oxide 99 layer thickness within the ps level. The initial oxide layer thickness, h, is estimated using a linear growth relationship, as shown in Equation 4.1: ℎ=¡¢u£ Equation 4.1 where C is related to the temperature dependent reaction rates. For the experimentally applied oxygen partial pressure of 10-5 Pa, it takes minutes to observe the oxide layer up to 10~100 nm thickness. Therefore ¡ = 10 A(cid:15)/(10(cid:22)⁄ ¢(cid:159)×100 ¥). With the same C, to form ~1nm oxide layer within 10~100 ps level in MD simulation, it requires PO2 = h/Ct = 1 nm/(10 nm×100 ps) × (10-5 Pa×100 s) = 10-6 Pa×1012 s = 106 Pa = 10 atm. Therefore, a much higher oxygen pressure was used. The oxidation rate under this condition is estimated by oxidizing pure Li slab in MD simulation. The number of oxygen atoms reacted with Li on the surface is 4~5 O atoms per 10 ps per nm2 under P = 20 atm. After adding oxygen molecules to the simulation cell the first time, the Li will be oxidized to Li oxide and consume oxygen molecules. To maintain a relative constant oxygen partial pressure, a special method is implemented in our simulation to maintain the constant oxygen partial pressure. First, we calculated the desired number of oxygen molecules needed in the simulation cell at different temperature and pressure. Then the oxygen molecules are filled into the simulation with the random position. After 4 ps of simulation, part of the oxygen molecule is reacted with surface Li. So we remove all unreacted oxygen molecule based on their charge and fill in new oxygen molecule to correct partial pressure again. In our simulation, the surface oxidation rate is about 0.5 ps-1 nm2 at 20 atm. 100 To compare with the morphology of oxidized Li surface with oxidized Al surface, a similar method is implemented to oxidize Al (0 0 1) slab for 250 ps. The simulation temperature is at scaled room temperature for Al. Details are explained in Chapter 4.4.1. 4.3.2. Extrusion of Li nanowire Figure 4.3 Schematic for slab model in MD simulation. The growth of Li nanowire is simulated by applying stress in oxygen environment. The oxidation rate is controlled by applying different temperature. To simulate the Li nanowire growth, extrusion of Li is applied to the slab model. As shown in Figure 4.3. the Li atoms in an 8 nm × 8 nm × 20 nm Li slab were extruded from a 2 nm diameter circle (dash region) on the surface, as the rest of the surface atoms and 1nm thick bottom (blue region) atoms are fixed. The fixed atoms on the top surface are subject to a constant velocity of - 0.2 nm/ps along normal of slab surface direction. In every 4ps, the velocity is applied for 2ps then 101 followed by 2 ps relaxation with surface position fixed. The atoms in 2 nm diameter circle as well as the bulk region (purple region) are subject to full relaxation for all simulation time. Other ways of simulating extrusion are also tested. For example, we applied a constant force on the surface atoms. But due to the diffusivity of Li atoms, the fully relaxed bulk Li atoms will penetrate the top surface when stress is large enough. We also tested the method to apply the stress by scaling the whole simulation cell along z-direction. But this will lead to continuously increasing of oxygen partial pressure during the extrusion. So, we used the extrusion method as introduced above. 4.3.3. Different conditions in MD simulation Corresponding to the different conditions used in the experiment, the oxidation rate and growth rate are carefully controlled in MD simulation. Three different conditions are designed in this study. For condition 1, the simulation temperature is 100 K with no additional oxygen. This temperature corresponds to 40% of Li melting temperature predicted in current ReaxFF forcefield (0.4Tm). At condition 1, the extrusion rate is ~ 243 Li atoms per 10 ps in nanowire, estimated from the growing nanoparticle from 25~35ps. For condition 2, the simulation temperature is lowered to 10 K with additional oxygen molecules of 20 atm partial pressure. The temperature is around 5% of Li melting temperature for a low oxidation rate. For condition 3, the oxygen partial pressure is set to 20 atm and the simulation temperature is 40% of Li melting temperature. Under this oxidation condition 3, the Li grew out were almost completely oxidized into oxide, reaching a Li:O ratio of 3:1. It is worth noting that, in condition 2, the Li extrusion rate is also slower than that in vacuum due to the oxide passivation. As a combined result, the Li growth rate ~5 times faster than the oxygen 102 atom incorporation rate, thus a thin surface oxide layer was formed. Indeed, the Li/O ratio in the nanowire in condition 2 is 5:1, since the center and the bottom of the NW is not completely oxidized. We have compared the following two cases: in case one, O2 pressure is applied for 10 ps until extrusion occurs. In the other case, the O2 pressure is applied 10 ps after the Li is extruded. Both cases show similar trends. 4.4. Results and discussion 4.4.1. Validation of ReaxFF reactive forcefield The accuracy of forcefield based molecular dynamics depends on the formulation and the corresponding parameters of the forcefield. The ReaxFF reactive forcefield is fitted to the Li-O system in previous works. [302] Due to the limitation of ReaxFF formulation, not all properties can be accurately predicted in the MD simulation. However, as long as the properties related to our research question are reasonably accurate, we can still gain useful insight from MD simulations. Therefore, the validation of forcefield is critical for the MD simulation. We calculated the melting point of Li using the current forcefield to be 250 K. It is lower than experimental value of 453.7 K because the melting point is not specifically included in the training data of fitting the forcefield. Although the absolute temperature is not accurate, we can still compare relative temperature and control the oxidation rate by changing the temperature in our MD simulation. [61] The melting point of Al is 670 K. So the absolute temperature of 500 K is used in the simulation, which has similar scale ratio as Li. 103 4.4.2. MD simulation of the oxidation process on Li and Al During the oxidation process (NVT dynamics under oxygen partial pressure of 10 atm and 150 K) the oxides gradually formed on the surface of Li slab. The top view and the side view of structures at different time step are shown in Figure 4.4. As a comparison, the Al morphology evolution is shown in Figure 4.5. The Al oxidation in MD simulation showed initially the oxides nucleate and grow non-uniformly, then form a continuous film, after which the oxide grows layer-by-layer and the thickness increases, which agrees well with the previous MD simulation [303]. A passivating oxide layer was formed with 1.5 nm thick after 250 ps of oxidation. This is also in consistency with experimental observations. The side view shows that the oxide layer is flat after oxidization. However, the Li surface showed a different pattern during oxidation. At the early stage (before 30 ps), the oxygen combines with Li atoms at the uniformly distributed location as shown in Figure 4.4. Then the oxides gradually aggregate to specific locations and become isolated clusters on the surface. Eventually, instead of forming a smooth passivating layer, a big crack is formed in the oxide layer. More cracks can be initiated from the location with thinner oxide layers under longer oxidation time. The side view shows the oxide has a non-uniform thickness and the Li-surface is not smooth anymore due to the deformation of the oxide layer. 104 Figure 4.4 Surface morphology evolution of Li slab under oxidation at different time step. For each structure, the top view and the side view is showed. The red atoms indicate oxygen and purple atoms indicate Li. Only Li and O in the oxides are shown for a clear view. 105 Figure 4.5 Surface morphology evolution of Al slab under oxidation at different time step. For each structure, the top view and the side view is showed. The red atoms indicate oxygen and grey atoms indicate Al. Only Al and O in the oxides are shown for a clear view. 106 The number of oxygen atoms reacted with metal in the two slabs is shown in Figure 4.6. At the early stage before 30 ps, the two slabs showed similar growth rate. After the Al slab surface being covered by oxides, the oxidation rate slows down following √§., indicating a diffusion controlled growth model after a passivation layer is formed. But the oxidation rate of Li maintains almost constant slope until 90 ps, indicating Li2O is non-passivating. Figure 4.6 Number of oxygen atoms reacted with metal during the oxidation process. The Al slab shows more passivating effect than Li. The oxidation simulation shows that Li2O is not an effective passivation layer compared to, for example, Al2O3, it shrinks in volume, cracks and exposes new pristine Li surface for oxidation. Cracking of oxide layer has also been observed at micron level experimentally. This interesting 107 observation is consistent with the empirical criteria for passivation oxide layer, known as the PB- ratio (PBR). The PBR was proposed by N.B. Pilling and R.E. Bedworth as an empirical relation for mechanical stability of simple surface metal oxides. It assumes stress generation in the surface oxide layer can cause crack formation, spallation of the surface oxide due to volume mismatch. Thus, PBR is defined as the molar volume ratio of the oxide and the metal. When PBR < 1, the coating layer may easily break; when PBR > 2, coating delaminates; when 1 < PBR < 2, the coating is passivating and provides a protecting effect against further surface oxidation. The density of Al oxides ranges from 3.5 to 4.0 g/cm3, which is 30% higher than the density of Al (2.7 g/cm3). So the PBR ratio for aluminum oxide is 1.3, (1<1.3<2). So a thin and dense passivating layer can form to prevent further oxidation. But the density of Li oxide is 2.0 g/cm3 while the density of Li is 0.5 g/cm3. After oxidation, the volume will shrink ~50%. So the PBR ratio is 0.6, less than 1. As expected, the lithium oxide layer will aggregate and lead to cracks on the surface. The opened cracks will expose fresh Li to oxygen. This explains the higher oxidation rate of Li comparing with Al. After 110 ps, the oxidation rate of Li also slows down because the Li oxide layer covered more and more surface area. Eventually, the Li oxides will also cover and passivate the Li surface. The Li oxides formed on the surface is curved in the end as shown in the side view in Figure 4.4. This is because, after the formation of thick Li oxides, the region covered by the oxide stopped reacting with oxygen. Li is only consumed from the cracks. Driven by the large Li oxides/Li interface energy, the newly formed oxides in the crack naturally diffused to the thick Li oxides region. So eventually, the previously covered region has a higher height along the surface normal direction. For Al slab, since the formation of Al oxides does not expose fresh Al, the whole surface can be oxidized at the same rate and result in a flat oxide layer. 108 4.4.3. Li growth morphology in oxygen environment The Li nanowire growth is simulated with different oxygen partial pressure and temperature. The cross-section of the final structures is shown in Figure 4.7. Li atoms are purple and oxygen atoms are yellow. For Figure 4.7(a) in condition 1, no extra oxygen is added and the simulation temperature is 40% of Li melting temperature. For Figure 4.7(b) in condition 2, oxygen with 20 atm partial pressure is added under the simulation temperature of 10K. For Figure 4.7(c) in condition 3, oxygen with 20 atm partial pressure is added under the simulation temperature of 40% Li melting temperature. The different growth mode of Li nanowire can be identified using the aspect ratio of the primitive cluster formed by Li atoms pushed out from the top surface. For condition 1, the Li formed a cluster with a height of 4.0 nm and diameter of 6.5 nm with aspect ratio less than 1 which indicates 3D growth mode. For condition 3, the height and diameter of Li cluster is 0.8 nm and 1.8 nm which as shows 3D growth mode. For condition 2, the height and diameter of Li cluster is 4.0 nm and 3.4 nm with aspect ratio larger than 1. It shows similar 1D growth mode in analogy to experimental observations. Also, we observed that a Li-core/oxide-shell structure is formed on the surface during the charging process as shown in Figure 4.7 (b). In the cross-section, most oxygen atoms aggregate on the outside layer of Li atoms which indicate the formation of Li oxides shell. These MD observations are consistent with experimental results discussed above. Specifically, in condition 2, the oxygen is consumed from the freshly exposed root of nanowire during the extrusion. This is due to the passivating effect of Li oxides on the outside of nanowire and will be further discussed in the next chapter. 109 Figure 4.7 MD simulation snapshots of Li grown under (a) UHV, (b) oxygen atmosphere, and (c) oxygen atmosphere at an increased temperature to reach a higher oxidation rate after 50 ps for parts a and c and 55 ps for part b, respectively The proposed mechanism for the effects of surface oxidation on the Li deposit morphology is depicted in Figure 4.8. We propose that it is surface oxidation of growing Li nanostructures and, more specifically, the interplay between the Li current densities and oxidation rate that are responsible for drastic changes in morphology of Li deposits: from 3D particles under the inert (Figure 4.7 a and Figure 4.8 a) and highly oxidizing ambient (Figure 4.7 c and Figure 4.8 c) to the formation of quasi-1D submicron thin nanowires when plating takes place at ≈ 6 × 10-6 Pa of oxygen partial pressure (Figure 4.7 b and Figure 4.8 b). The formation of Li-core/oxide-shell structure around growing Li structures is known from previous studies[301]. AES measurements[151] also show brightness reversal of corresponding SEM images due to higher secondary electron yield of Li oxide and its surface charging. The development of Li-oxide shells around the growing metallic nanowires is further evidenced via SEM imaging[151] of the same 110 anode region before and after Li discharge. The SEM images show bright, solid shells with a thickness of few tens of nanometers remaining after the metallic Li core is intercalated back into LiCoO2 during battery discharge. Figure 4.8 Schematics showing models of the Li nucleation and growth for the images presented in Figure 4.7. In the presence of O2, a lithium-oxide sheath is developed. Schematics are drawn not to scale. For an intermediate oxygen pressure range ≈ 6 × 10-6 Pa, the O atomic flux rate is ≈ 5 × 1013 (s·cm2)-1, which is 100 times smaller than the Li current density ¤. Therefore, the Li nucleation at the anode surface proceeds initially similarly to UHV conditions. The difference however is that the growing Li surface become immediately oxidized. The thickness h of the newly formed oxide depends of the Li ion diffusion D in the oxide (Li ions/vacancies diffuse significantly faster than oxygen ones), oxygen partial pressure PO2 and exposure time t [304]: 111 ℎ=j2'“¢u£a(cid:23) § Equation 4.2 here N1 is the number of oxygen atoms in Li oxide unit cell and K is the Henry’s law constant. The presence of a thin oxide sheath prevents the surface from rapid bulk oxidation and inhibits Li surface diffusion[305,306] due to the high activation energy. Therefore, the volume expansion during charging can only occur only at the footprint of the Li particles through the out-of-plane elongation resulting in lifting of particle out-of-plane (Figure 4.7 b). As the Li nanowire grow, newly emerged Li surface at the wire footprint gets immediately oxidized and encapsulated, thus promoting a one-dimensional growth (Figure 4.7 b) similar to the vapor-liquid-solid growth mechanism of metal oxide nanowires or the oxygen assisted growth of Si nanowires. Similarly, a thin oxide layer is formed also at the surface of Li saturated carbon anode. This eventually prevents unimpeded Li diffusion towards the free anode surface and explains the observed lower number density of Li nuclei. Since more Li will be accumulated inside the C-anode, a developing compressive stress is channeling Li to existing growth spots. Finally, an additional driving force for 1D growth develops due to progressive Li nanowire surface oxidation. Due to the higher density of Li2O compare to metallic Li the surface oxide layer shrinks in volume and compressive stress is developed inside NW. The stress is released via the directed flow of the metallic Li core towards the NW footprint. This 1D Li growth proceeds until partial oxygen pressure reaches the level where a sufficiently thick oxide layer is formed around the Li nuclei and at the surface of the lithiated carbon. At this pressure, the thickness of the oxide near the Li particle footprint reaches the critical thickness 112 ℎ(cid:136)f«≥›fi(cid:158)fl(cid:16)(cid:176) Equation 4.3 where it will not fracture upon the given compressive stress in the C-anode, here ›fi is the shear stress at the Li/Li2O interface, (cid:16)(cid:176) is the oxide fracture strength, (cid:158)fl is the distance of the Li growth front to the stress-free end of the oxide layer.[307] Note, this mechanically robust oxide shell is formed around each cluster including regions close to the root of the Li nanostructures (Figure 4.7 c), pinning Li deposits to the C anode surface and inhibiting any further growth. Since the potential barrier for particle nucleation is proportional to γ3, where γ is the surface energy of Li, oxidation of Li leads to surface energy increase of Li-O interface, resulting in potential barrier growth. Therefore, the Li flux J contributes mostly to nucleation of new Li nanostructures which are more energetically favorable rather than to the growth of the pre-existing ones, in a manner analogous to metal electroplating from solution in the presence of brightening agent such as saccharin.[302] One can still use scaling Equation 4.2 to explain qualitatively an increase of nucleation density at elevated oxygen partial pressures. An oxide formation at the C-anode surface implies slower Li delivery rate J and as a result, a longer time for oxide to form on the surface. Therefore, higher N is observed under high oxidant partial pressure. 4.4.4. Competition between oxidation rates and growth rates The 1D nanowire only formed under the balanced charge and oxidation rates. Although Figure 4.7 only showed one rate and three oxygen partial pressures. One can expect the effect of the charge rate. A slow charge rate means more time for oxide to form on the surface. Therefore one may see 3D growth under lower oxygen partial pressure and lower charging rate. On the other hand, a faster 113 charging rate means less time for oxide layer to form, therefore one may see 1D wire formation under the high oxygen partial pressure but faster rate. Therefore, we predict Li morphology changes from microparticles to Nanowires and back to microparticles under a combined effect of charging rates and oxygen partial pressure, in Table 4.1. Table 4.1 The expected Li growth morphology at different charging rate and oxygen pressure. The growth mechanism is observed based on Li nanowire growing on the carbon anode surface. But we expect it can be generalized to other materials with similar structures such as Sn whisker growth on Sn coating [308]. The impact of internal stress imposed by misorientation between grains and external oxidation effect may have similar competition relationship that can be described by the model in this study. 114 4.5. Conclusion Molecule dynamics model is built to model a thin film solid state battery system and provided new insights into the mechanisms of early stage Li plating and the surprising effects of oxidant on the Li nanostructure growth. The experimental study shows that under an inert (UVH) or relatively high (in excess of ≈ 10-5 Pa) oxygen partial pressure low rate (≤1 mA/cm2) Li plating proceeds via nucleation and growth of 3D particles. However, at intermediate O2 partial pressure of ≈ 10-6 Pa, the oxide sheath formed around the growing Li nanostructures promotes the Li deposit growth in the form of 1D nanowires. This transition is largely controlled by the thickness of the oxide layer formed on the Li surface and therefore depends strongly on the ratio between charging (J) and oxide thickness growth ((cid:139)–(cid:139)(cid:134)) rates. We anticipate that such conditions favoring 1D growth can exist in a wide range of J·((cid:139)–(cid:139)(cid:134))(cid:22)(cid:23)ratios. Very fast charging, however, leads to Li segregation inside the volume of the anode, compressive stress development, cracking/delamination followed by Li metal creeping at the interfaces and large whiskers growth independent of ambient conditions. 115 CHAPTER 5. Improvement of charge transfer scheme in ReaxFF 5.1. Summary The forcefield based molecular dynamics method (MD) can bring significant insights into understanding deformation mechanisms. Specifically, the qEq method based ReaxFF reactive forcefield has been successfully used in metal and ceramic systems to predict atomic charges. However, it also has limitations in predicting atomic charges, and consequently leads to different simulated fracture behavior. In qEq method, the instantaneous charge transfer and long-range charge exchange is allowed. This introduces extra energy even the bond breaks in the ionic system, which will lead to an over-ductile problem in fracture simulation. To expose this issue, we demonstrated the over-ductile problem in the qEq scheme by simulating tensile of a Li2O slab model. It clearly showed the Li oxides with brittle nature undergo ductile fracture during the deformation. To fix this problem, we proposed to replace the qEq method with ACKS2 method. The ACKS2 parameters are obtained from high-level quantum calculations. The impact of charge calculation method on the fracture behavior of Li2O slab is analyzed. 5.2. Forcefield and its limitation on fracture simulation Forcefield-based Molecular dynamics (MD) simulation has been widely used to explore the dislocation interactions, grain boundary sliding, and fracture behavior of nano-structured materials [309–312]. Despite their ability to accurately describe the equilibrium structure and unveil nanoscale deformation mechanisms, traditional forcefield methods fail to determine the ductile and brittle fracture behavior in many cases. For example, Kang implemented five different forcefields for the tensile test of Si and Ge nanowire. Among them, three potentials predicted 116 ductile fracture mode while the other two potentials predicted brittle fracture mode. Also, the predicted ideal tensile strength for Si nanowire can vary from 10.2 GPa to 26.3 GPa depending on the choice of forcefield [312]. The observation of ductile or brittle fracture behavior varies with the geometry of the structures. Islam et al. observed the ductile fracture of Li0.8S using a nanowire structure without notch. The fracture is initiated via formation of voids predicted using ReaxFF reactive forcefield [313]. Brutzel et al. studied the fracture of silica in a notched bulk structure and observed brittle fracture using SW potential [314]. Even if the simulation is carried out by the same potential on the same material, the predicted results for toughness may be different based on the structure geometry. Yuan compared the tensile test of bulk silica with and without a notch in the structure both using BKS potential. The MD results showed that the elongation to break can be 0.23 without notch or 0.11 with a notch [315]. This over-ductile behavior in MD simulations is caused by various simulation conditions such as the size of the system, the geometry of the structure, and the choice of forcefields. If the size is too small, it is difficult to nucleate any defects such as voids and dislocations. So the material will act as a perfect single crystal with higher yield strength, fracture strength and more brittleness than the experimental behavior. Even when the system size is large enough, due to the limitation of simulation time, the strain rate is much higher than the typical tensile test which may lead to large unphysical elongation to break in a perfect structure because the defects do not have enough time to nucleate and propagate. Additional notch on the structure can act as pre-exist defects to facilitate the nucleation of initial cracks. So, a notch should be added to the structure for fracture simulation. 117 Describing chemical bond formation and breaking properly remain challenging. In the formulation of some forcefields, important interactions might be ignored due to simplicity. Fracture behavior is a bond-breaking process. For example, in metallic systems, the embedded atom potential (EAM) ignored the directional bonds due to the free electron ocean assumption. The interaction with far neighbors is approximated by interactions with uniformly distributed local electron density. Therefore, EAM potential intrinsically cannot describe to highly anisotropic materials well. The Modified EAM potential (MEAM) added angle contributions and showed improvements. For ionic materials, many simple forcefields with fixed charges cannot lack the ability to describe the charge transfer process accurately and dynamically. Notably, during large deformation and fracture, when two atoms are far apart, they should become neutral species rather than charged. It is anticipated that the residue charge on the atoms after bond breaking will lead to extra Coulomb interactions which overestimate the toughness. Several methods have been used to allocate the atomic charge dynamically. Sanderson [316] stated the principle of electronegativity equalization that electrons will flow until electronegativities on all atoms are equalized. Mortier et al. [317,318] developed the electronegativity equalization method (EEM) that minimizes the charge energy as a quadratic function of charge distribution. After Mortier, various modifications have been added to the EEM schemes. Rappe and Goddard [319] proposed the qEq method based on EEM as shown in Equation 5.1: V=V‡·(cid:181)¶•+.l‚0‚6 4 „ (cid:128)(cid:19)p¤(|>0,6|) 0,6 +.Χ0‚0 0 +.»02‚0(cid:24) 0 Equation 5.1 118 where the first term is short-ranged interaction, the second term is screened Coulomb interaction, Χ0 is the electronegativity and »0 is the atomic hardness. The qEq method is widely used in forcefields that involve charge transfer including ReaxFF [133] and COMB [132]. The qEq method based ReaxFF reactive forcefield has been used in many studies to understand the deformation process [320–322]. Ostadhossein et al. developed ReaxFF potential for Li-Si-O- Al system and simulated the Li diffusion in silica [323]. Kim et al. implemented this potential to investigate the fracture of the coating during the lithiation process of Si-core/Al2O3-shell and Si- core/SiO2-shell nanoclusters. The compositional change, stress distribution and fracture behavior under different lithiation rate are simulated [143]. However, it is realized that the charge prediction in the current ReaxFF method is not always accurate. It is often difficult to determine if a phenomenon is due to real physics or the limitation of forcefields. In ionic systems, the qEq method suffers from unphysical long-range interactions and instantaneous charge transfer. So it is implicitly assumed the charge can transfer from one atom to another instantaneously. In other words, all materials are conductive based on this treatment. Also, in the energy minimization process, no constrains are imposed on the amount of charge that can be transferred. This means even with large distance apart, two atoms may still exchange and carry fractional charge when they supposed to be neutral species. The non-neutrality leads to additional long-range Coulomb interactions which eventually result in the over-ductile of ionic materials, which typically are brittle. To fix this problem, some ad hoc fixes are proposed such as adding artificial constraints on atomic charges (CHARMM) [324]. A different approach called atom-atom charge transfer (AACT) method proposed by Chelli [325] introduced a dummy variable, split charge, to share charge across 119 a bond. The split charge determines the amount of charge transferred from one atom to the other. Later, Nistor et al. [326] proposed split-charge equilibration (SQE) method which combines the EEM and AACT method. The formalism is shown as: V=V‡·(cid:181)¶•+.l‚0‚6 4 „ (cid:128)(cid:19)p¤(|>0,6|) 0,6 +.Χ0‚0 0 +.»02‚0(cid:24) 0 + . »0,62 …0,6(cid:24) 0,6,670 Equation 5.2 where qH,(cid:190) is the split charge associated with the bond between atom i and atom j. The SQE and its improved version can accurately capture some charge transfer phenomenon such as static electricity transfer and voltage change in battery system [327]. The atom-condensed Kohn-Sham DFT approximated to second order (ACKS2) method is another generation of SQE and qEq method. It is derived by expanding the Legendre transform of the Kohn-Sham kinetic energy to second order in the atomic populations and a new set of dual atomic variables, the relative atomic Kohn-Sham potentials [153]. It is proved that ACKS2 method will give a better description of atomic charge as well as a more accurate energy prediction [249]. As illustrated in Figure 5.1, the charge on hydrogen in the HF molecule is predicted using three different methods as a function of dissociation distance. The result of ACKS2 is consistent with the high-level quantum mechanics calculation (CAS-SCF/6-311++G**) at the intermediate distance as well as long distance. The qEq can only predict the correct charge at the equilibrium distance as indicated by the vertical line in the figure. At long dissociation distance, the EEM failed to predict the neutrality of atoms. 120 Figure 5.1 Numerical results of charge on hydrogen atom for the dissociation of hydrogen fluoride. [153] In this chapter, the research goal is to explore the impact of the charge transfer method on the fracture behavior in simple ionic materials. The fracture behavior of a notched Li2O slab under tension will be simulated and the crack propagation patterns will be tracked to determine the ductile or brittle fracture characteristics. First, the qEq charge allocation method along with the other parameters developed for the ReaxFF will be tested to explore the simulated fracture behavior. Then, the more advanced charge transfer method, ACKS2 method with the new forcefield parameters fitted to high accuracy quantum calculations using a genetic algorithm (GA), will be used to simulate the fracture behavior. It is anticipated that the ACKS2 may mitigate the over-ductile problem faced by ReaxFF in the ionic system since the ACKS2 can eliminate the residual charge problem during bond breaking. This chapter is conducted in collaboration with Dr. Aktulga’s group [154,155]. 121 5.3. Computational approach The qEq based ReaxFF MD simulations are carried out using LAMMPS package[251,252]. The qEq-based ReaxFF parameter used in this study is developed by Alireza Ostadhossein et al. [323]. The ACSK2 based ReaxFF MD simulations are carried out using PuReMD code developed by Aktulga et al. The parameters are newly fitted based on DFT calculated energetics for various Li and Li2O crystal structures and the atomic charge during bond breaking based on high-level quantum chemistry calculations conducted by Illias Magoulas in Dr. Piotr’s group. The training data includes the equation of states of Li2O crystal in BCC, FCC, simple cubie, diamond and HCP structures; the dissociation of Li-Li and Li-O-Li molecules; the shearing and tensile of Li2O and LiO crystal. DFT calculations for additional training data are conducted by Michael Swift from MSCE group in MSU. The ACKS2 parameter fitting is conducted by Kurt O’Hearn from Prof. Metin Aktulga in CSE department in MSU. Figure 5.2 The (a) as build structure of Li2O. Li atoms are in purple. Oxygen atoms are in red. 122 The Li2O unit cell structure with Fm3m symmetry is obtained from reference [328]. A crystal slab model with 12672 Li and 6336 O atoms is created from the unit cell. The cell length is 18.44 Å× 100 Å×82.98 Å. To create the slab surface, 40 Å vacuum is added along the y-direction. Finally, a 30-degree notch with length of 25 Å along surface normal direction is added to the surface of the slab as shown in Figure 5.2. Minimization of the slab model is done by applying NVT dynamics at 1K for 20000 steps with 0.25 fs time step for both methods. The relaxed structures using the qEq method and ACKS2 method are shown in Figure 5.3. The slab model is stretched along the z-direction. The cell length along z-direction is increased by 0.5% of its original length followed by 1 ps of NVT dynamics for relaxation. After stretching, the atomic positions are mapped to the stretched cell such that their fractional coordinates remain the same. This is equivalent to a strain rate of 5×10(cid:192)/s. The tensile tests are simulated at 1K and 200K to compare the impact of temperature. 5.4. Results and discussion 5.4.1. Relaxation of the slab model Figure 5.3 Relaxed Li2O slab structure using (a) qEq and (b) ACKS2 method. Li atoms are in purple. Oxygen atoms are in red. 123 Figure 5.4 Radial distribution function of two Li2O crystal structures and two relaxed slab structures As shown in Figure 5.3(a), after 20000 steps of relaxation at 1 K using the qEq method, the crystal structure has transformed into the amorphous structure. For ACKS2 relaxed slab model in Figure 5.3(b), a more ordered structure is maintained comparing with the qEq method. This improvement is due to the adjustment of forcefield parameters in our fitting process. The challenging part of this unexpected phase transition might be due to the fact that two crystalline structures of Li2O exist. The stable structure is cubic with Fm-3m symmetry and the structure with PNMA symmetry is only less than 0.1 eV higher than the cubic Li2O. Such small energy difference between the two 124 crystalline structures is generally hard to distinguish. Never the less, the radial distribution function of Li2O crystal and slab structures shown in Figure 5.4 showed that apart from the shift of the first peak from 2.0 Å to 1.7 Å, both slab structures showed a peak at 2.4 Å which is similar with crystal structures. The ACKS2 relaxed slab also has a broad peak at 2.0 Å which appeared in the crystal structures. 5.4.2. Deformation of the notched Li2O slab model Figure 5.5 Deformation of Li2O slab model using qEq method under (a) 1K and (b) 200K at different strains. Deformation with lower temperature showed more ductility during fracture. 125 The deformation process of the slab model at 1K and 200K using the qEq method is shown in Figure 5.5. At 200 K, the fracture front showed a blunt tip which indicates a ductile fracture. This is contradicting with the brittleness nature of ceramic Li2O. As we discussed above, this over- ductile phenomenon is a systemic error in the current qEq method due to the long-range interaction caused by incorrect prediction of atomic charge. When the deformation is conducted under 1K which will significantly increase the brittleness of any material, the crack tip showed brittle fracture features as shown in Figure 5.5 (a). This means although not all predicted properties are correct, the forcefield method can still capture general physics behaviors during the fracture process. The stress-strain curve of these two deformation tests shown in Figure 5.6 supported the previous conclusion. The 1K tensile test (blue curve) showed higher fracture strength and sudden drop of stress around 20% strain, which clearly shows brittle fracture characteristics. While the 200 K tensile test (red curve) showed much higher elongation to break with gradual decreasing of stress value after yield which indicates ductile fracture behavior. 126 Figure 5.6 The stress-strain curve of Li2O slab deformation process using the qEq method at 1K and 200K 5.5. Future work Next, we will keep on working on improving the fitting of forcefield parameters and conduct more detailed analysis on the tensile test of Li2O slab model with the same procedure using ACKS2 method. The stress-strain curve will be plotted and compared with qEq results. From the shape of the crack tip and stress-strain curve, the brittleness of Li2O slab can be estimated. By comparing the predicted brittleness using qEq and ACKS2 method, we can see the effect of eliminating long- range interaction on the previous over-ductile problem. 127 CHAPTER 6. Summary and outlook In this thesis, we implemented the atomic simulation techniques including DFT and forcefield based MD methods to explore the deformation mechanism of complex nano-structured materials. Also, we demonstrated one possible solution for the over-ductile problem in the current charge transfer scheme. To understand the deformation mechanism of nacre under dynamic loading, the fracture energy, generalized stacking fault energy curve, and twinning tendency of (110)[1 10] slip system is calculated using DFT method. Because of the atomic relaxation along (110) normal direction in nacre structure, the unstable stacking fault energy falls below fracture energy of mineral aragonite which avoids the fracture of nacre under impact. The twinning tendency result predicts the dynamic impact will favor formation of twinning under dynamic loading condition rather than full dislocation, which explains the experimental observations. The unstable/stable stacking fault energy of nacre shows the highest value among a wide range of metals and ceramics. The low stable stacking fault energy is due to its large number of multi-neighbor shared ionic bonds that provided similar chemical environment before and after sliding. The high unstable stacking fault energy is due to the oxygen-oxygen repulsion and the unique relaxation process during the sliding. Overall, the deformation mechanism corresponding to nacre’s nanoscale structure works together with other macroscale deformation mechanisms in synergistically provide the extraordinary toughening and strengthening effect under impact loading. To understand the impact of formation history and chemical-mechanical coupled effect from the environment on the fracture behavior of aluminum bifilm, a ReaxFF reactive forcefield based MD 128 procedure is established to simulate its fracture process and mechanical properties. Our results show that during the aging process, the fracture energy of bifilm increased from 0.43 J/m2 to 0.53 J/m2 for young and old bifilm. The clean “dry” interface of bifilm can “heal” itself to 50% to 75% of the bonding density comparing with bulk amorphous Al2O3. This effect leads to fracture at “wetted” interface, Al/Al2O3 interface, of bifilm in vacuum for both young and old bifilm. However, when the environmental factor is considered, the “dry” interface of bifilm can be contaminated by hydrogen due to the rapid solubility change at a different temperature. With 30% hydroxyl group coverage, the bifilm tend to fracture at the Al2O3/Al2O3 interface. The fracture strength of young and old bifilm are 1.1 GPa and 1.55 GPa. With surface contamination, the fracture strength drops to 1.2 GPa. Therefore, our simulation results indicate that the crack initiation position is influenced by surface contamination more than the phase change from the aging process. So effective removal of hydrogen gas may lead to a strengthening of bifilm interface and better mechanical performance of casting product. Also, we converted the MD simulated fracture strength to macroscale cohesive zone model parameters for FEM simulation. To understand the growth mechanism of Li dendrite under different oxygen partial pressure, we implemented ReaxFF reactive forcefield based MD simulation to model the competition between oxidation rate and dendrite growth rate. The results show that the MD method can simulate the early stage of Li dendrite growth on carbon anode. When oxygen is in deficient (10-7 Pa), the Li will diffuse into all directions on the anode surface which lead to the formation of large clusters. When oxygen partial pressure is high enough, a thick Li oxide layer will form on the surface and suppress the nucleation of Li dendrite due to the high modulus of Li oxide. This will lead also lead to diffusion of Li along all directions and form small clusters on the surface. When the oxidation rate and growth rate are compatible, a core-shell structure is formed that provided support and 129 confinement effect on the in-plane direction, which leads to preferred growth on the plane normal direction. The MD simulation validated the proposed growth mechanism and demonstrated the importance of trace amount of oxygen on the morphology evolution of the lithiation process in all- solid-state Li battery. Finally, to improve the current charge transfer scheme used in forcefield based MD simulation, we first expose the over-ductile issue in the current method using tensile test of Li2O slab. It shows unphysical ductile fracture behavior. The parameters for ACKS2 method are fitted from highly accurate training data obtained by high-level quantum chemistry calculations. The fracture of Li2O slab with the same procedure using the ACKS2 method shows it can alleviate the over-ductile problem in the fracture process due to the more accurate atomic charge prediction. Overall, we demonstrate the capability of various atomic simulation techniques on predicting mechanical properties, fracture behaviors, and most importantly, the deformation mechanisms of nano-structured materials. The atomic simulation techniques show superior ability over experimental techniques on observing and measuring the chemical and mechanical process simultaneously for nano-structured material applications. To correctly capture the mechanical process and environmental effect in our model, we carefully designed simulation protocols based on experimental conditions. These protocols are successfully implemented to bring insights into the deformation mechanisms of metal/ceramics composite materials with nanoscale substructures. It opens up many opportunities to further explore the impact of chemical reaction on the deformation process of other nano-structured materials. Also, we realize the current simulation tools are not perfect. With the rapid development of new mathematics tools and deeper physical understanding of materials, new models should be 130 continuously developed to improve the accuracy and efficiency of atomic simulation methods. For example, the recently blooming of machine learning method in computer vision and national language processing enabled many new tools to unleash the power of big data and modern supercomputers. Specifically, the rapid development of deep learning techniques has achieved better than human-level ability on pattern recognition tasks with the help of rational designed neural net structures and a large amount of computational power. This technique can be used in high-throughput calculations for composition and structure screening to find materials with desired properties. On the other hand, generative adversarial model (GAN) have the ability to learn underlying data structure and probability distribution from real data and generate similar data from the same distribution [330]. The conditional GAN can take a “condition” of each training data as input and learn the relation between input condition and distribution of training data [331]. This method can be used as a surrogate model of the real physics-based model but with enormous amount of reduced evaluation time. These physics informed deep learning model could be used to replace one or more steps of traditional high computational cost models such as solving for turbulence models in computational fluid dynamics (CFD) [332] and atomic charge distribution in density functional theory [333]. 131 BIBLIOGRAPHY 132 BIBLIOGRAPHY [1] W.F. Smith, J. Hashemi, Foundation of materials science and engineering, 4th ed., McGraw-Hill, New York, 2006. [2] C. Mercer, W.O. Soboyejo, Hall-petch relationships in gamma titanium aluminides, Scr. Mater. 35 (1996) 17–22. doi:10.1016/1359-6462(96)00097-8. [3] Y.-H. Kim, A. Inoue, T. Masumoto, Increase in Mechanical Strength of Al-Y-Ni Amorphous Alloys by Dispersion of Nanoscale fcc-Al Particles, Mater. Trans. JIM. 32 (1991) 331–338. doi:10.2320/matertrans1989.32.331. [4] [5] T. Namazu, Y. Isono, T. Tanaka, Evaluation of size effect on mechanical properties of single crystal silicon by nanoscale bending test using AFM, J. Microelectromechanical Syst. 9 (2000) 450–459. doi:10.1109/84.896765. E.W. Qin, L. Lu, N.R. Tao, J. Tan, K. Lu, Enhanced fracture toughness and strength in bulk nanocrystalline Cu with nanoscale twin bundles, Acta Mater. 57 (2009) 6215–6225. doi:10.1016/J.ACTAMAT.2009.08.048. [6] H. Zhou, S. Qu, The effect of nanoscale twin boundaries on fracture toughness in doi:10.1088/0957- 035706. (2010) 21 nanocrystalline Ni, Nanotechnology. 4484/21/3/035706. [7] J.E. Allen, E.R. Hemesath, D.E. Perea, J.L. Lensch-Falk, Z.Y. Li, F. Yin, M.H. Gass, P. Wang, A.L. Bleloch, R.E. Palmer, L.J. Lauhon, High-resolution detection of Au catalyst atoms in Si nanowires, Nat. Nanotechnol. 3 (2008) 168–173. doi:10.1038/nnano.2008.5. [8] Y. Cui, C.M. Lieber, M.G. Bawendi, Functional nanoscale electronic devices assembled 851–3. Science. building blocks., (2001) 291 using doi:10.1126/science.291.5505.851. nanowire silicon [9] C. Lee, X. Wei, J.W. Kysar, J. Hone, Measurement of the elastic properties and intrinsic strength of monolayer graphene., Science. 321 (2008) 385–8. doi:10.1126/science.1157996. [10] A. Fasolino, J.H. Los, M.I. Katsnelson, Intrinsic ripples in graphene, Nat. Mater. 6 (2007) 858–861. doi:10.1038/nmat2011. [11] A. Castellanos-Gomez, M. Poot, G.A. Steele, H.S.J. van der Zant, N. Agraït, G. Rubio- Bollinger, Elastic Properties of Freely Suspended MoS 2 Nanosheets, Adv. Mater. 24 (2012) 772–775. doi:10.1002/adma.201103965. [12] S. Bertolazzi, J. Brivio, A. Kis, Stretching and Breaking of Ultrathin MoS 2, ACS Nano. 5 (2011) 9703–9709. doi:10.1021/nn203879f. 133 [13] I. Jo, M.T. Pettes, J. Kim, K. Watanabe, T. Taniguchi, Z. Yao, L. Shi, Thermal Conductivity and Phonon Transport in Suspended Few-Layer Hexagonal Boron Nitride, Nano Lett. 13 (2013) 550–554. doi:10.1021/nl304060g. [14] Y. Li, Y. Rao, K.F. Mak, Y. You, S. Wang, C.R. Dean, T.F. Heinz, Probing Symmetry Properties of Few-Layer MoS 2 and h-BN by Optical Second-Harmonic Generation, Nano Lett. 13 (2013) 3329–3333. doi:10.1021/nl401561r. [15] B. Anasori, M.R. Lukatskaya, Y. Gogotsi, 2D metal carbides and nitrides (MXenes) for energy storage, Nat. Rev. Mater. 2 (2017) 16098. doi:10.1038/natrevmats.2016.98. [16] A. Gupta, T. Sakthivel, S. Seal, Recent development in 2D materials beyond graphene, Prog. Mater. Sci. 73 (2015) 44–126. doi:10.1016/J.PMATSCI.2015.02.002. [17] R. Mas-Ballesté, C. Gómez-Navarro, J. Gómez-Herrero, F. Zamora, 2D materials: to graphene and beyond, Nanoscale. 3 (2011) 20–30. doi:10.1039/C0NR00323A. [18] W. Han, Perspectives for spintronics in 2D materials, APL Mater. 4 (2016) 032401. doi:10.1063/1.4941712. [19] H. Li, Z.H. Xu, X. Li, Multiscale hierarchical assembly strategy and mechanical prowess in 409–416. J. Struct. Biol. carica), (2013) 184 conch doi:10.1016/j.jsb.2013.10.011. (Busycon shells [20] M. Holzinger, A. Le Goff, S. Cosnier, Nanomaterials for biosensing applications: a review, Front. Chem. 2 (2014) 63. doi:10.3389/fchem.2014.00063. [21] S. Guo, E. Wang, Noble metal nanomaterials: Controllable synthesis and application in fuel 240–264. analytical Today. (2011) Nano 6 cells sensors, doi:10.1016/J.NANTOD.2011.04.007. and [22] M.U. Niemann, S.S. Srinivasan, A.R. Phani, A. Kumar, D.Y. Goswami, E.K. Stefanakos, Nanomaterials for Hydrogen Storage Applications: A Review, J. Nanomater. 2008 (2008) 1–9. doi:10.1155/2008/950967. [23] M.D. Uchic, D.M. Dimiduk, J.N. Florando, W.D. Nix, Sample dimensions influence strength and crystal plasticity, Science. 305 (2004) 986–989. doi:10.1126/science.1098993. [24] Y. Kondo, K. Takayanagi, Synthesis and characterization of helical multi-shell gold nanowires, Science. 289 (2000) 606–8. doi:10.1126/SCIENCE.289.5479.606. [25] Y. Kondo, K. Takayanagi, Gold Nanobridge Stabilized by Surface Structure, Phys. Rev. Lett. 79 (1997) 3455–3458. doi:10.1103/PhysRevLett.79.3455. [26] H. Ohnishi, Y. Kondo, K. Takayanagi, Quantized conductance through individual rows of suspended gold atoms, Nature. 395 (1998) 780–783. doi:10.1038/27399. [27] D. François, A. Pineau, A. Zaoui, Mechanical Behaviour of Materials, Springer Netherlands, 134 Dordrecht, 2013. doi:10.1007/978-94-007-4930-6. [28] V. Richter, M.. Ruthendorf, On hardness and toughness of ultrafine and nanocrystalline hard materials, Int. J. Refract. Met. Hard Mater. 17 (1999) 141–152. doi:10.1016/S0263- 4368(99)00003-7. [29] Y. Wang, M. Chen, F. Zhou, E. Ma, High tensile ductility in a nanostructured metal, Nature. 419 (2002) 912–915. doi:10.1038/nature01133. [30] Y.S. Sato, M. Urata, H. Kokawa, K. Ikeda, Hall–Petch relationship in friction stir welds of equal channel angular-pressed aluminium alloys, Mater. Sci. Eng. A. 354 (2003) 298–305. doi:10.1016/S0921-5093(03)00008-X. [31] A.H. Chokshi, A. Rosen, J. Karch, H. Gleiter, On the validity of the hall-petch relationship in nanocrystalline materials, Scr. Metall. 23 (1989) 1679–1683. doi:10.1016/0036- 9748(89)90342-6. [32] A.M. El-Sherik, U. Erb, G. Palumbo, K.T. Aust, Deviations from hall-petch behaviour in (1992) 1185–1188. as-prepared nanocrystalline nickel, Scr. Metall. Mater. 27 doi:10.1016/0956-716X(92)90596-7. [33] N. Hansen, Hall–Petch relation and boundary strengthening, Scr. Mater. 51 (2004) 801–806. doi:10.1016/J.SCRIPTAMAT.2004.06.002. [34] P.G. Sanders, J.A. Eastman, J.R. Weertman, Elastic and tensile behavior of nanocrystalline copper and palladium, Acta Mater. 45 (1997) 4019–4025. doi:10.1016/S1359- 6454(97)00092-X. [35] A. Giga, Y. Kimoto, Y. Takigawa, K. Higashi, Demonstration of an inverse Hall–Petch relationship in electrodeposited nanocrystalline Ni–W alloys through tensile testing, Scr. Mater. 55 (2006) 143–146. doi:10.1016/J.SCRIPTAMAT.2006.03.047. [36] P.G. Sanders, J.A. Eastman, J.R. Weertman, Elastic and tensile behavior of nanocrystalline copper and palladium, Acta Mater. 45 (1997) 4019–4025. doi:10.1016/S1359- 6454(97)00092-X. [37] V.Y. Gertsman, M. Hoffmann, H. Gleiter, R. Birringer, The study of grain size dependence of yield stress of copper for a wide grain size range, Acta Metall. Mater. 42 (1994) 3539– 3544. doi:10.1016/0956-7151(94)90486-3. [38] M. Zhao, J.C. Li, Q. Jiang, Hall–Petch relationship in nanometer size range, J. Alloys Compd. 361 (2003) 160–164. doi:10.1016/S0925-8388(03)00415-8. [39] A.G. Evans, J.P. Hirth, Deformation of nanoscale cermets, Scr. Metall. Mater. 26 (1992) 1675–1680. doi:10.1016/0956-716X(92)90532-J. [40] C. Suryanarayana, D. Mukhopadhyay, S.N. Patankar, F.H. Froes, Grain size effects in nanocrystalline materials, J. Mater. Res. 7 (1992) 2114–2118. doi:10.1557/JMR.1992.2114. 135 [41] D.A. Konstantinidis, E.C. Aifantis, On the `anomalous’ hardness of nanocrystalline (1998) 1111–1118. doi:10.1016/S0965- materials, Nanostructured Mater. 10 9773(98)00145-7. [42] Y. Qi, P.E. Krajewski, Molecular dynamics simulations of grain boundary sliding: The effect of stress and boundary misorientation, Acta Mater. 55 (2007) 1555–1563. doi:10.1016/j.actamat.2006.10.016. [43] O.D. Sherby, J. Wadsworth, Superplasticity—Recent advances and future directions, Prog. Mater. Sci. 33 (1989) 169–221. doi:10.1016/0079-6425(89)90004-2. [44] K. Higashi, High strain rate superplasticity in Japan, Mater. Sci. Technol. 16 (2000) 1320– http://www.scopus.com/inward/record.url?eid=2-s2.0- 1329. 0000291661&partnerID=40&md5=b17d07ddce506d4944585f1f8dcd3daf. [45] S.X. McFadden, R.S. Mishra, R.Z. Valiev, A.P. Zhilyaev, A.K. Mukherjee, Low- temperature superplasticity in nanostructured nickel and metal alloys, Nature. 398 (1999) 684–686. doi:10.1038/19486. [46] R.S. Mishra, S.X. McFadden, R.Z. Valiev, A.K. Mukherjee, Deformation mechanisms and (1999) 37–40. in nanocrystalline materials, JOM. 51 tensile doi:10.1007/s11837-999-0010-1. superplasticity [47] R.S. Mishra, T.R. Bieler, A.K. Mukherjee, Mechanism of high strain rate superplasticity in aluminium alloy composites, Acta Mater. 45 (1997) 561–568. doi:10.1016/S1359- 6454(96)00194-2. [48] L. Zhang, C. Lu, K. Tieu, A review on atomistic simulation of grain boundary behaviors in 180–191. cubic metals, Comput. Mater. Sci. (2016) face-centered doi:10.1016/J.COMMATSCI.2016.03.021. 118 [49] X. Li, Z.H. Xu, R. Wang, In situ observation of nanograin rotation and deformation in nacre, Nano Lett. 6 (2006) 2301–2304. doi:10.1021/nl061775u. [50] D.F. Hou, G.S. Zhou, M. Zheng, Conch shell structure and its effect on mechanical behaviors, Biomaterials. 25 (2004) 751–756. doi:10.1016/S0142-9612(03)00555-6. [51] Z. Huang, Z. Pan, H. Li, Q. Wei, X. Li, Hidden energy dissipation mechanism in nacre, J. Mater. Res. 29 (2014) 1573–1578. doi:10.1557/jmr.2014.179. [52] J. Sun, B. Bhushan, Hierarchical structure and mechanical properties of nacre: a review, RSC Adv. 2 (2012) 7617. doi:10.1039/c2ra20218b. [53] K.S. Katti, D.R. Katti, S.M. Pradhan, A. Bhosle, Platelet interlocks are the key to toughness and strength in nacre, J. Mater. Res. 20 (2005) 1097–1100. doi:10.1557/JMR.2005.0171. [54] B.L. Smith, T.E. Schaffer, M. Viani, J.B. Thompson, N. a Frederick, J. Kindt, a M. Belcher, G.D. Stucky, D.E. Morse, P.K. Hansma, Molecular mechanistic origin of the toughness of 136 natural adhesives, fibres and composites, Nature. 399 (1999) 761–763. doi:10.1038/21607. [55] I.R. Kramer, L.J. Demer, Effects of environment on mechanical properties of metals, Prog. Mater. Sci. 9 (1961) 133–199. doi:10.1016/0079-6425(61)90020-2. [56] J.K. Chang, E.M. Taleff, P.E. Krajewski, J.R. Ciulik, Effects of atmosphere in filament formation on a superplastically deformed aluminum-magnesium alloy, Scr. Mater. 60 (2009) 459–462. doi:10.1016/j.scriptamat.2008.11.031. [57] W.H. Krueger, S.R. Pollack, The initial oxidation of aluminum thin films at room temperature, Surf. Sci. 30 (1972) 263–279. doi:10.1016/0039-6028(72)90002-7. [58] F. Jona, Preparation and properties of clean surfaces of aluminum, J. Phys. Chem. Solids. 28 (1967) 2155–2160. doi:10.1016/0022-3697(67)90239-9. [59] F.G. Sen, A.T. Alpas, A.C.T. Van Duin, Y. Qi, Oxidation-assisted ductility of aluminium nanowires, Nat. Commun. 5 (2014) 3959. doi:10.1038/ncomms4959. [60] Y. Yang, A. Kushima, W. Han, H. Xin, J. Li, Liquid-Like, Self-Healing Aluminum Oxide during Deformation at Room Temperature, Nano Lett. 18 (2018) 2492–2497. doi:10.1021/acs.nanolett.8b00068. [61] J. Liu, Q. Wang, Y. Qi, Atomistic simulation of the formation and fracture of oxide bifilms 673–682. (2019) Mater. cast in Acta doi:10.1016/J.ACTAMAT.2018.11.008. aluminum, 164 [62] K. Xu, Electrolytes and Interphases in Li-Ion Batteries and Beyond, Chem. Rev. 114 (2014) 11503–11618. doi:10.1021/cr500003w. [63] X.-B. Cheng, R. Zhang, C.-Z. Zhao, F. Wei, J.-G. Zhang, Q. Zhang, A Review of Solid Electrolyte Interphases on Lithium Metal Anode, Adv. Sci. 3 (2016) 1500213. doi:10.1002/advs.201500213. [64] Z.H. Xu, Y. Yang, Z. Huang, X. Li, Elastic modulus of biopolymer matrix in nacre measured using coupled atomic force microscopy bending and inverse finite element techniques, Mater. Sci. Eng. C. 31 (2011) 1852–1856. doi:10.1016/j.msec.2011.08.023. [65] E.C.C.M. Silva, L. Tong, S. Yip, K.J. Van Vliet, Size Effects on the Stiffness of Silica Nanowires, Small. 2 (2006) 239–243. doi:10.1002/smll.200500311. [66] A. Heidelberg, L.T. Ngo, B. Wu, M.A. Phillips, S. Sharma, T.I. Kamins, J.E. Sader, J.J. Boland, A generalized description of the elastic properties of nanowires, Nano Lett. 6 (2006) 1101–1106. doi:10.1021/nl060028u. [67] B. Wu, A. Heidelberg, J.J. Boland, Mechanical properties of ultrahigh-strength gold nanowires, Nat. Mater. 4 (2005) 525–529. doi:10.1038/nmat1403. [68] D.M. Kochmann, From atoms to the macroscale: the nonlocal quasicontinuum method, 137 (2019). http://www.mm.ethz.ch/research_QC.html. [69] S. Iijima, High‐Resolution Electron Microscopy of Crystal Lattice of Titanium‐Niobium Oxide, J. Appl. Phys. 42 (1971) 5891–5893. doi:10.1063/1.1660042. [70] J.M. Cowley, S. Iijima, Electron Microscope Image Contrast for Thin Crystal, Zeitschrift Für Naturforsch. A. 27 (1972) 445–451. doi:10.1515/zna-1972-0312. [71] Z. Horita, D.J. Smith, M. Furukawa, M. Nemoto, R.Z. Valiev, T.G. Langdon, An investigation of grain boundaries in submicrometer-grained Al-Mg solid solution alloys using high-resolution electron microscopy, J. Mater. Res. 11 (1996) 1880–1890. doi:10.1557/JMR.1996.0239. [72] T.J. Balk, K.J. Hemker, High resolution transmission electron microscopy of dislocation core dissociations in gold and iridium, Philos. Mag. A. 81 (2001) 1507–1531. doi:10.1080/01418610108214360. [73] J.M. Gibson, M.L. McDonald, F.C. Unterwald, Direct imaging of a novel silicon surface reconstruction, Phys. Rev. Lett. 55 (1985) 1765–1767. doi:10.1103/PhysRevLett.55.1765. [74] H. Hashimoto, T. Naiki, T. Eto, K. Fujiwara, High Temperature Gas Reaction Specimen Chamber for an Electron Microscope, Jpn. J. Appl. Phys. 7 (1968) 946–952. doi:10.1143/JJAP.7.946. [75] H. Yoshida, S. Takeda, T. Uchiyama, H. Kohno, Y. Homma, Atomic-Scale In-situ Observation of Carbon Nanotube Growth from Solid State Iron Carbide Nanoparticles, Nano Lett. 8 (2008) 2082–2086. doi:10.1021/nl080452q. [76] S. Helveg, C. López-Cartes, J. Sehested, P.L. Hansen, B.S. Clausen, J.R. Rostrup-Nielsen, F. Abild-Pedersen, J.K. Nørskov, Atomic-scale imaging of carbon nanofibre growth, Nature. 427 (2004) 426–429. doi:10.1038/nature02278. [77] M. Li, D.-G. Xie, E. Ma, J. Li, X.-X. Zhang, Z.-W. Shan, Effect of hydrogen on the integrity of aluminium–oxide interface at elevated temperatures, Nat. Commun. 8 (2017) 14564. doi:10.1038/ncomms14564. [78] D. Xie, S. Li, M. Li, Z. Wang, P. Gumbsch, J. Sun, E. Ma, J. Li, Z. Shan, Hydrogenated (2016) 13341. in aluminium, Nat. Commun. 7 vacancies doi:10.1038/ncomms13341. lock dislocations [79] M.T. Reetz, W. Helbig, S.A. Quaiser, U. Stimming, N. Breuer, R. Vogel, Visualization of Surfactants on Nanostructured Palladium Clusters by a Combination of STM and High- Resolution TEM., Science. 267 (1995) 367–9. doi:10.1126/science.267.5196.367. [80] F. Wang, H.-C. Yu, M.-H. Chen, L. Wu, N. Pereira, K. Thornton, A. Van der Ven, Y. Zhu, G.G. Amatucci, J. Graetz, Tracking lithium transport and electrochemical reactions in nanoparticles, Nat. Commun. 3 (2012) 1201. doi:10.1038/ncomms2185. 138 [81] M. Gu, L.R. Parent, B.L. Mehdi, R.R. Unocic, M.T. McDowell, R.L. Sacci, W. Xu, J.G. Connell, P. Xu, P. Abellan, X. Chen, Y. Zhang, D.E. Perea, J.E. Evans, L.J. Lauhon, J.-G. Zhang, J. Liu, N.D. Browning, Y. Cui, I. Arslan, C.-M. Wang, Demonstration of an Electrochemical Liquid Cell for Operando Transmission Electron Microscopy Observation of the Lithiation/Delithiation Behavior of Si Nanowire Battery Anodes, Nano Lett. 13 (2013) 6106–6112. doi:10.1021/nl403402q. [82] X.H. Liu, J.W. Wang, S. Huang, F. Fan, X. Huang, Y. Liu, S. Krylyuk, J. Yoo, S.A. Dayeh, A. V. Davydov, S.X. Mao, S.T. Picraux, S. Zhang, J. Li, T. Zhu, J.Y. Huang, In situ atomic- scale imaging of electrochemical lithiation in silicon, Nat. Nanotechnol. 7 (2012) 749–756. doi:10.1038/nnano.2012.170. [83] J.Y. Huang, L. Zhong, C.M. Wang, J.P. Sullivan, W. Xu, L.Q. Zhang, S.X. Mao, N.S. Hudak, X.H. Liu, A. Subramanian, H. Fan, L. Qi, A. Kushima, J. Li, In situ observation of the electrochemical lithiation of a single SnO‐ nanowire electrode., Science. 330 (2010) 1515–20. doi:10.1126/science.1195628. [84] Y. Li, Y. Li, A. Pei, K. Yan, Y. Sun, C.-L. Wu, L.-M. Joubert, R. Chin, A.L. Koh, Y. Yu, J. Perrino, B. Butz, S. Chu, Y. Cui, Atomic structure of sensitive battery materials and interfaces revealed by cryo-electron microscopy., Science. 358 (2017) 506–510. doi:10.1126/science.aam6014. [85] J.A. Pople, R. Krishnan, H.B. Schlegel, J.S. Binkley, Derivative studies in hartree-fock and 225–241. J. Quantum Chem. (2009) Int. 16 møller-plesset doi:10.1002/qua.560160825. theories, [86] C. David Sherrill, H.F. Schaefer, The Configuration Interaction Method: Advances in (1999) 143–269. Highly Correlated Approaches, Adv. Quantum Chem. 34 doi:10.1016/S0065-3276(08)60532-8. [87] R.J. Bartlett, M. Musiał, Coupled-cluster theory in quantum chemistry, Rev. Mod. Phys. 79 (2007) 291–352. doi:10.1103/RevModPhys.79.291. [88] S. Grimme, Improved second-order Møller–Plesset perturbation theory by separate scaling of parallel- and antiparallel-spin pair correlation energies, J. Chem. Phys. 118 (2003) 9095– 9102. doi:10.1063/1.1569242. [89] W. Kohn, L.J. Sham, Self-Consistent Equations Including Exchange and Correlation Effects, Phys. Rev. 140 (1965) A1133–A1138. doi:10.1103/PhysRev.140.A1133. [90] P. Hohenberg, W. Kohn, Inhomogeneous Electron Gas, Phys. Rev. 136 (1964) B864–B871. doi:10.1103/PhysRev.136.B864. [91] M. P. Allen, D.J. Tildesley, Computer simulation of liquid metals, 2nd ed., Oxford University Press, New York, 2017. [92] E.B. Tadmor, R.E. Miller, Modeling materials: Continuum, atomistic and multiscale 2011. Cambridge, techniques, Cambridge University Press, 139 doi:10.1017/CBO9781139003582. [93] J.L. Bouvard, D.K. Ward, D. Hossain, S. Nouranian, E.B. Marin, M.F. Horstemeyer, Review of Hierarchical Multiscale Modeling to Describe the Mechanical Behavior of Amorphous Polymers, J. Eng. Mater. Technol. 131 (2009) 041206. doi:10.1115/1.3183779. [94] H. Van Swygenhoven, P.M. Derlet, a G. Frøseth, Stacking fault energies and slip in nanocrystalline metals., Nat. Mater. 3 (2004) 399–403. doi:10.1038/nmat1136. [95] G. Lu, N. Kioussis, V. V. Bulatov, E. Kaxiras, Generalized-stacking-fault energy surface and dislocation properties of aluminum, Phys. Rev. B. 62 (2000) 3099–3108. doi:10.1103/PhysRevB.62.3099. [96] V. Yamakov, D. Wolf, S.R. Phillpot, H. Gleiter, Deformation twinning in nanocrystalline (2002) 5005–5020. Al by molecular dynamics simulation, Acta Mater. 50 doi:10.1016/S1359-6454(02)00318-X. [97] V. Yamakov, D. Wolf, S.R. Phillpot, H. Gleiter, Dislocation – dislocation and dislocation – twin reactions in nanocrystalline Al by molecular dynamics simulation, Acta Mater. 51 (2003) 4135–4147. doi:10.1016/S1359-6454(03)00232-5. [98] E.B. Tadmor, S. Hai, A Peierls criterion for the onset of deformation twinning at a crack tip, J. Mech. Phys. Solids. 51 (2003) 765–793. doi:10.1016/S0022-5096(03)00005-X. [99] G. Lu, N. Kioussis, V. V Bulatov, E. Kaxiras, Generalized stacking fault energy surfaces (1999) 25. aluminum, Phys. Rev. B. 62 and dislocation properties of doi:10.1103/PhysRevB.62.3099. [100] M.J. Cawkwell, D. Nguyen-Manh, C. Woodward, D.G. Pettifor, V. Vitek, Origin of brittle cleavage in iridium., Science. 309 (2005) 1059–62. doi:10.1126/science.1114704. [101] C. Woodward ∥, S.I. Rao, Ab-initio simulation of ( a /2)⟨110] screw dislocations in γ-TiAl, Philos. Mag. 84 (2004) 401–413. doi:10.1080/14786430310001611626. [102] C. Woodward, S.I. Rao, Flexible Ab Initio Boundary Conditions: Simulating Isolated (2002) 216402. in bcc Mo and Ta, Phys. Rev. Lett. 88 Dislocations doi:10.1103/PhysRevLett.88.216402. [103] C. Woodward, D.R. Trinkle, L.G. Hector, D.L. Olmsted, Prediction of Dislocation Cores in Aluminum from Density Functional Theory, Phys. Rev. Lett. 100 (2008) 045507. doi:10.1103/PhysRevLett.100.045507. [104] L.A. Zepeda-Ruiz, A. Stukowski, T. Oppelstrup, V. V. Bulatov, Probing the limits of metal (2017) 492–495. plasticity with molecular dynamics simulations, Nature. 550 doi:10.1038/nature23472. [105] M. Levitt, The birth of computational structural biology, Nat. Struct. Biol. 8 (2001) 392– 393. https://www.nature.com/articles/nsb0501_392 (accessed January 20, 2019). 140 [106] M.S. Daw, M.I. Baskes, Embedded-atom method: Derivation and application to impurities, in metals, Phys. Rev. B. 29 (1984) 6443–6453. surfaces, and other defects doi:10.1103/PhysRevB.29.6443. [107] J.D. Gale, Empirical potential derivation for ionic materials, Philos. Mag. B. 73 (1996) 3– 19. doi:10.1080/13642819608239107. [108] S.J. Weiner, P.A. Kollman, D.T. Nguyen, D.A. Case, An all atom force field for simulations (1986) 230–252. J. Comput. Chem. 7 and nucleic acids, of proteins doi:10.1002/jcc.540070216. [109] S.J. Weiner, P.A. Kollman, D.A. Case, U.C. Singh, C. Ghio, G. Alagona, S. Profeta, P. Weiner, A new force field for molecular mechanical simulation of nucleic acids and proteins, J. Am. Chem. Soc. 106 (1984) 765–784. doi:10.1021/ja00315a051. [110] D.Q. McDonald, W.C. Still, AMBER torsional parameters for the peptide backbone, Tetrahedron Lett. 33 (1992) 7743–7746. doi:10.1016/0040-4039(93)88034-G. [111] A.K. Rappe, C.J. Casewit, K.S. Colwell, W.A. Goddard, W.M. Skiff, UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations, J. Am. Chem. Soc. 114 (1992) 10024–10035. doi:10.1021/ja00051a040. [112] C.J. Casewit, K.S. Colwell, A.K. Rappe, Application of a universal force field to main group compounds, J. Am. Chem. Soc. 114 (1992) 10046–10053. doi:10.1021/ja00051a042. [113] C.J. Casewit, K.S. Colwell, A.K. Rappe, Application of a universal force field to organic molecules, J. Am. Chem. Soc. 114 (1992) 10035–10046. doi:10.1021/ja00051a041. [114] A.K. Rappe, K.S. Colwell, C.J. Casewit, Application of a universal force field to metal complexes, Inorg. Chem. 32 (1993) 3438–3450. doi:10.1021/ic00068a012. [115] S.L. Mayo, B.D. Olafson, W.A. Goddard, DREIDING: a generic force field for molecular simulations, J. Phys. Chem. 94 (1990) 8897–8909. doi:10.1021/j100389a010. [116] D. Farkas, W.A. Curtin, Plastic deformation mechanisms in nanocrystalline columnar grain structures, Mater. Sci. Eng. A. 412 (2005) 316–322. doi:10.1016/J.MSEA.2005.09.043. [117] D.E. Spearot, K.I. Jacob, D.L. McDowell, Nucleation of dislocations from [0 0 1] bicrystal 3579–3589. (2005) Mater. Acta interfaces doi:10.1016/J.ACTAMAT.2005.04.012. aluminum, in 53 [118] A. Luque, M. Ghazisaeidi, W.A. Curtin, Deformation modes in magnesium (0 0 0 1) and $(0\,1\,\bar{1}\,1)$ single crystals: simulations versus experiments, Model. Simul. Mater. Sci. Eng. 21 (2013) 045010. doi:10.1088/0965-0393/21/4/045010. [119] T. KITAMURA, K. YASHIRO, R. OHTANI, Atomic Simulation on Deformation and Fracture of Nano-Single Crystal of Nickel in Tension., JSME Int. J. Ser. A. 40 (1997) 430– 435. doi:10.1299/jsmea.40.430. 141 [120] V. Yamakov, D. Wolf, S.R. Phillpot, A.K. Mukherjee, H. Gleiter, Dislocation processes in the deformation of nanocrystalline aluminium by molecular-dynamics simulation., Nat. Mater. 1 (2002) 45–8. doi:10.1038/nmat700. [121] D. Rodney, Molecular dynamics simulation of screw dislocations interacting with interstitial frank loops in a model FCC crystal, Acta Mater. 52 (2004) 607–614. doi:10.1016/j.actamat.2003.09.044. [122] V. Yamakov, D. Wolf, S.R. Phillpot, H. Gleiter, Grain-boundary diffusion creep in nanocrystalline palladium by molecular-dynamics simulation, Acta Mater. 50 (2002) 61– 73. doi:10.1016/S1359-6454(01)00329-9. [123] J.Q. Broughton, G.H. Gilmer, Thermodynamic criteria for grain-boundary melting: A 2692–2695. (1986) Phys. Rev. Lett. molecular-dynamics doi:10.1103/PhysRevLett.56.2692. study, 56 [124] J. Lao, M. Naghdi Tam, D. Pinisetty, N. Gupta, Molecular Dynamics Simulation of FCC Metallic Nanowires: A Review, JOM. 65 (2013) 175–184. doi:10.1007/s11837-012-0465- 3. [125] J. Diao, K. Gall, M. Dunn, Surface stress driven reorientation of gold nanowires, Phys. Rev. B. 70 (2004) 075413. doi:10.1103/PhysRevB.70.075413. [126] J. Diao, K. Gall, M. L. Dunn, Atomistic simulation of the structure and elastic properties of 1935–1962. J. Mech. Solids. (2004) Phys. 52 gold doi:10.1016/J.JMPS.2004.03.009. nanowires, [127] J. Lao, D. Moldovan, Surface stress induced structural transformations and pseudoelastic effects in palladium nanowires, Appl. Phys. Lett. 93 (2008) 093108. doi:10.1063/1.2976434. [128] J. Diao, K. Gall, M.L. Dunn, Surface-stress-induced phase transformation in metal nanowires, Nat. Mater. 2 (2003) 656–660. doi:10.1038/nmat977. [129] V. Yamakov, D. Wolf, S.R. Phillpot, A.K. Mukherjee, H. Gleiter, Deformation-mechanism map for nanocrystalline metals by molecular-dynamics simulation, Nat. Mater. 3 (2004) 43–47. doi:10.1038/nmat1035. [130] T.S. Bush, J.D. Gale, C.R. a. Catlow, P.D. Battle, Self-consistent interatomic potentials for the simulation of binary and ternary oxides, J. Mater. Chem. 4 (1994) 831. doi:10.1039/jm9940400831. [131] M.F. Russo, A.C.T. Van Duin, Atomistic-scale simulations of chemical reactions: Bridging from quantum chemistry to engineering, Nucl. Instruments Methods Phys. Res. Sect. B Beam Interact. with Mater. Atoms. 269 (2011) 1549–1554. doi:10.1016/j.nimb.2010.12.053. [132] T. Liang, T.-R. Shan, Y.-T. Cheng, B.D. Devine, M. Noordhoek, Y. Li, Z. Lu, S.R. Phillpot, S.B. Sinnott, Classical atomistic simulations of surfaces and heterogeneous interfaces with the charge-optimized many body (COMB) potentials, Mater. Sci. Eng. R Reports. 74 (2013) 142 255–279. doi:10.1016/J.MSER.2013.07.001. [133] A.C.T. Van Duin, S. Dasgupta, F. Lorant, W.A. Goddard, ReaxFF: A reactive force field for hydrocarbons, J. Phys. Chem. A. 105 (2001) 9396–9409. doi:10.1021/jp004368u. [134] F.G. Sen, Y. Qi, A.C.T. Van Duin, A.T. Alpas, Oxidation induced softening in Al nanowires, Appl. Phys. Lett. 102 (2013) 051912. doi:10.1063/1.4790181. [135] B. Narayanan, A.C.T. van Duin, B.B. Kappes, I.E. Reimanis, C. V. Ciobanu, Reactive force field for lithium-aluminum silicates with applications to eucryptite phases, Model. Simul. Mater. Sci. Eng. 20 (2012) 015002. doi:10.1088/0965-0393/20/1/015002. [136] Kevin D. Nielson, Adri C. T. van Duin, Jonas Oxgaard, and Wei-Qiao Deng, W.A.G. III*, Development of the ReaxFF Reactive Force Field for Describing Transition Metal Catalyzed Reactions, with Application to the Initial Stages of the Catalytic Formation of Carbon Nanotubes, (2004). doi:10.1021/JP046244D. [137] Q. Zhang, T. Çağln, A. Van Duin, W.A. Goddard, Y. Qi, L.G. Hector, Adhesion and nonwetting-wetting transition in the Al/α−Al2O3interface, Phys. Rev. B - Condens. Matter Mater. Phys. 69 (2004) 045423. doi:10.1103/PhysRevB.69.045423. [138] X.G. Wang, J.R. Smith, A. Evans, Fundamental Influence of C on Adhesion of the 286102. (2002) Rev. Lett. Al2O3/Al Phys. doi:10.1103/PhysRevLett.89.286102. Interface, 89 [139] A. Sazgar, M.R. Movahhedy, M. Mahnama, S. Sohrabpour, A molecular dynamics study of bond strength and interface conditions in the Al/Al2O3 metal–ceramic composites, Comput. Mater. Sci. 109 (2015) 200–208. doi:10.1016/J.COMMATSCI.2015.07.024. [140] M. Yoshio, H. Wang, K. Fukuda, T. Umeno, N. Dimov, Z. Ogumi, Carbon-Coated Si as a Lithium-Ion Battery Anode Material, J. Electrochem. Soc. 149 (2002) A1598. doi:10.1149/1.1518988. [141] S.-H. Ng, J. Wang, D. Wexler, K. Konstantinov, Z.-P. Guo, H.-K. Liu, Highly Reversible Lithium Storage in Spheroidal Carbon-Coated Silicon Nanocomposites as Anodes for Lithium-Ion 7050–7053. doi:10.1002/ange.200601676. Batteries, Chemie. Angew. (2006) 118 [142] S. Sim, P. Oh, S. Park, J. Cho, Critical thickness of SiO2 coating layer on core@Shell bulk@nanowire Si anode materials for Li-ion batteries, Adv. Mater. 25 (2013) 4498–4503. doi:10.1002/adma.201301454. [143] S.-Y. Kim, A. Ostadhossein, A.C.T. van Duin, X. Xiao, H. Gao, Y. Qi, Self-generated concentration and modulus gradient coating design to protect Si nano-wire electrodes during lithiation, Phys. Chem. Chem. Phys. 18 (2016) 3706–3715. doi:10.1039/C5CP07219K. [144] X. Wu, D. Ma, P. Eisenlohr, D. Raabe, H.-O. Fabritius, From insect scales to sensor design: modelling the mechanochromic properties of bicontinuous cubic structures, Bioinspir. 143 Biomim. 11 (2016) 045001. doi:10.1088/1748-3190/11/4/045001. [145] D.D. Tjahjanto, P. Eisenlohr, F. Roters, Multiscale deep drawing analysis of dual-phase steels using grain cluster-based RGC scheme, Model. Simul. Mater. Sci. Eng. 23 (2015) 045005. doi:10.1088/0965-0393/23/4/045005. [146] Z. Huang, H. Li, Z. Pan, Q. Wei, Y.J. Chao, X. Li, Uncovering high-strain rate protection mechanism in nacre, Sci. Rep. 1 (2011) 1–5. doi:10.1038/srep00148. [147] Q.G. Wang, C.J. Davidson, J.R. Griffiths, P.N. Crepeau, Oxide films, pores and the fatigue lives of cast aluminum alloys, Metall. Mater. Trans. B Process Metall. Mater. Process. Sci. 37 (2006) 887–895. doi:10.1007/BF02735010. [148] X. Dai, X. Yang, J. Campbell, J. Wood, Effects of runner system design on the mechanical strength of Al–7Si–Mg alloy castings, Mater. Sci. Eng. A. 354 (2003) 315–325. doi:10.1016/S0921-5093(03)00021-2. [149] X. Yang, X. Huang, X. Dai, J. Campbell, J. Tatler, Numerical modelling of entrainment of oxide film defects in filling of aluminium alloy castings, Int. J. Cast Met. Res. 17 (2004) 321–331. doi:10.1179/136404604225022748. [150] C.M. Pita, S.D. Felicelli, A fluid–structure interaction method for highly deformable solids, Comput. Struct. 88 (2010) 255–262. doi:10.1016/J.COMPSTRUC.2009.11.004. [151] A. Yulaev, V. Oleshko, P. Haney, J. Liu, Y. Qi, A.A. Talin, M.S. Leite, A. Kolmakov, From Microparticles to Nanowires and Back: Radical Transformations in Plated Li Metal Morphology Revealed via in Situ Scanning Electron Microscopy, Nano Lett. 18 (2018) 1644–1650. doi:10.1021/acs.nanolett.7b04518. [152] N. Legrand, B. Knosp, P. Desprez, F. Lapicque, S. Raël, Physical characterization of the charging process of a Li-ion battery and prediction of Li plating by electrochemical modelling, 208–216. doi:10.1016/J.JPOWSOUR.2013.06.130. Sources. (2014) J. Power 245 [153] T. Verstraelen, P.W. Ayers, V. Van Speybroeck, M. Waroquier, ACKS2: Atom-condensed Kohn-Sham DFT approximated to second order, J. Chem. Phys. 138 (2013) 074108. doi:10.1063/1.4791569. [154] H.M. Aktulga, J.C. Fogarty, S.A. Pandit, A.Y. Grama, Parallel reactive molecular dynamics: Numerical methods and algorithmic techniques, Parallel Comput. 38 (2012) 245–259. doi:10.1016/J.PARCO.2011.08.005. [155] H.M. Aktulga, C. Knight, P. Coffman, K.A. O’Hearn, T.-R. Shan, W. Jiang, Optimizing the Performance of Reactive Molecular Dynamics Simulations for Multi-Core Architectures, (2017). http://arxiv.org/abs/1706.07772 (accessed February 21, 2019). [156] A.P. Jackson, J.F. V. Vincent, R.M. Turner, The Mechanical Design of Nacre, in: Proc. R. Soc. B Biol. Sci., 1988: pp. 415–440. doi:10.1098/rspb.1988.0056. 144 [157] J.H.E. Cartwright, A.G. Checa, The dynamics of nacre self-assembly., J. R. Soc. Interface. 4 (2007) 491–504. doi:10.1098/rsif.2006.0188. [158] S. Kamat, X. Su, R. Ballarini, a H. Heuer, Structural basis for the fracture toughness of the shell of the conch Strombus gigas., Nature. 405 (2000) 1036–1040. doi:10.1038/35016535. [159] H. Imai, Y. Oaki, The Hierarchical Architecture of Nacre and its Mimetic Materials, Handb. Biominer. Biol. Asp. Struct. Form. 2 (2008) 89–107. doi:10.1002/9783527619443.ch29. [160] G. Zhang, X. Li, Uncovering aragonite nanoparticle self-assembly in nacre-A natural armor, Cryst. Growth Des. 12 (2012) 4306–4310. doi:10.1021/cg3010344. [161] T.E. Schäffer, C. Ionescu-Zanetti, R. Proksch, M. Fritz, D.A. Walters, N. Almqvist, C.M. Zaremba, A.M. Belcher, B.L. Smith, G.D. Stucky, D.E. Morse, P.K. Hansma, Does Abalone Nacre Form by Heteroepitaxial Nucleation or by Growth through Mineral Bridges?, Chem. Mater. 9 (1997) 1731–1740. doi:10.1021/cm960429i. [162] F. Song, X.H. Zhang, Y.L. Bai, Microstructure and Characteristics in the Organic Matrix Layers of Nacre, J. Mater. Res. 17 (2002) 1567–1570. doi:10.1557/JMR.2002.0233. [163] Q.L. Feng, F.Z. Cui, G. Pu, R.Z. Wang, H.D. Li, Crystal orientation, toughening mechanisms and a mimic of nacre, Mater. Sci. Eng. C. 11 (2000) 19–25. doi:10.1016/S0928-4931(00)00138-7. [164] J.D. Currey, Mechanical Properties of Mother of Pearl in Tension, Proc. R. Soc. London. Ser. B, Biol. Sci. 196 (1977) 443–463. doi:10.1098/rspb.1977.0050. [165] A.P. Jackson, J.F. V Vincent, R.M. Turner, Comparison of nacre with other ceramic composites, J. Mater. Sci. 25 (1990) 3173–3178. doi:10.1007/BF00587670. [166] R.Z. Wang, Z. Suo, A.G. Evans, N. Yao, I.A. Aksay, Deformation mechanisms in nacre, J. Mater. Res. 16 (2001) 2485–2493. doi:10.1557/JMR.2001.0340. [167] F. Barthelat, C.-M. Li, C. Comi, H.D. Espinosa, Mechanical properties of nacre constituents and their impact on mechanical performance, J. Mater. Res. 21 (2006) 1977–1986. doi:10.1557/jmr.2006.0239. [168] F. Barthelat, H. Tang, P.D. Zavattieri, C.M. Li, H.D. Espinosa, On the mechanics of mother- of-pearl: A key feature in the material hierarchical structure, J. Mech. Phys. Solids. 55 (2007) 306–337. doi:10.1016/j.jmps.2006.07.007. [169] K.S. Katti, D.R. Katti, Why is nacre so tough and strong?, Mater. Sci. Eng. C. 26 (2006) 1317–1324. doi:10.1016/j.msec.2005.08.013. [170] B. Ji, H. Gao, A study of fracture mechanisms in biological nano-composites via the virtual 96–103. bond model, Mater. (2004) internal doi:10.1016/j.msea.2003.08.121. Sci. Eng. A. 366 145 [171] R.Z. Wang, H.B. Wen, F.Z. Cui, H.B. Zhang, H.D. Li, Observations of damage morphologies in nacre during deformation and fracture, J. Mater. Sci. 30 (1995) 2299–2304. doi:10.1007/BF01184577. [172] D.R. Katti, K.S. Katti, Modeling microarchitecture and mechanical behavior of nacre using 3D finite element techniques. Part I. Elastic properties, J. Mater. Sci. 36 (2001) 1411–1417. doi:10.1023/A:1017528209162. [173] H.H. Fu, D.J. Benson, M.A. Meyers, Analytical and computational description of effect of grain size on yield stress of metals, Acta Mater. 49 (2001) 2567–2582. doi:10.1016/S1359- 6454(01)00062-3. [174] H. Mukai, K. Saruwatari, H. Nagasawa, T. Kogure, Aragonite twinning in gastropod nacre, J. Cryst. Growth. 312 (2010) 3014–3019. doi:10.1016/j.jcrysgro.2010.07.002. [175] Y.A. Shin, S. Yin, X. Li, S. Lee, S. Moon, J. Jeong, M. Kwon, S.J. Yoo, Y.-M. Kim, T. Zhang, H. Gao, S.H. Oh, Nanotwin-governed toughening mechanism in hierarchically structured biological materials., Nat. Commun. 7 (2016) 10772. doi:10.1038/ncomms10772. [176] S.N. Wang, X.H. Yan, R. Wang, D.H. Yu, X.X. Wang, A microstructural study of individual nacre tablet of Pinctada maxima, J. Struct. Biol. 183 (2013) 404–411. doi:10.1016/j.jsb.2013.07.013. [177] X. Li, W.C. Chang, Y.J. Chao, R. Wang, M. Chang, Nanoscale structural and mechanical characterization of a natural nanocomposite material: The shell of red abalone, Nano Lett. 4 (2004) 613–617. doi:10.1021/nl049962k. [178] T.Z. Blazynski, Materials at High Strain Rates, ELSEVIER APPLIED SCIENCE PUBLISHERS LTD, 1987. https://books.google.co.il/books?id=1rb7TRrrAM8C. [179] J.D. Campbell, W.G. Ferguson, The temperature and strain-rate dependence of the shear strength of mild steel, Philos. Mag. 21 (1970) 63–82. doi:10.1080/14786437008238397. [180] A.R. Dowling, J. Harding, J.D. Campbell, The dynamic punching of metals, J. Inst. Met. 98 (1970) 215–224. [181] Y.L. Bai, J. Harding, Fracture initiation in glass-reinforced plastics under impact loading, Mechanical Properties at High Rates of Strain, in: Inst. Phys. Conf., Oxford, 1984. [182] P.A.A. Back, J.D. Campbell, The Behaviour of a Reinforced Plastic Material Under Dynamic Compression, in: Proc. Conf. Prop. Mater. High Rates Strain, 1957: pp. 221–228. [183] V. Vítek, Intrinsic stacking faults in body-centred cubic crystals, Philos. Mag. 18 (1968) 773–786. doi:10.1080/14786436808227500. [184] D.J. Siegel, Generalized stacking fault energies, ductilities, and twinnabilities of Ni and selected Ni alloys, Appl. Phys. Lett. 87 (2005) 121901. doi:10.1063/1.2051793. 146 [185] X.L. Wu, Y. Qi, Y.T. Zhu, Partial-mediated slips in nanocrystalline Ni at high strain rate, Appl. Phys. Lett. 90 (2007) 22–24. doi:10.1063/1.2745250. [186] H. Jónsson, G. Mills, K.W. Jacobsen, Nudged elastic band method for finding minimum energy paths of transitions, World Scientific, 1998. doi:10.1142/9789812839664_0016. [187] G. Henkelman, H. Jónsson, Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points, J. Chem. Phys. 113 (2000) 9978–9985. doi:10.1063/1.1323224. [188] G. Henkelman, G. Jóhannesson, H. Jónsson, Theoretical Methods in Condensed Phase Chemistry, in: Theor. Methods Condens. Phase Chem. - Prog. Theor. Chem. Phys. - Vol. 5, Springer, 2002: pp. 269–302. doi:10.1007/0-306-46949-9. [189] D. Sheppard, R. Terrell, G. Henkelman, Optimization methods for finding minimum energy paths, J. Chem. Phys. 128 (2008) 134106. doi:10.1063/1.2841941. [190] DFT_method, No Title, (2016). [191] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations (1996) 11169–11186. set, Phys. Rev. B. 54 a plane-wave basis using doi:10.1103/PhysRevB.54.11169. [192] G. Kresse, J. Hafner, Ab liquid- metalamorphous- semiconductor transition in germanium, Phys. Rev. B. 49 (1994) 14251– 14269. doi:10.1103/PhysRevB.49.14251. initio molecular-dynamics simulation of the [193] G. Kresse, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B. 59 (1999) 1758–1775. doi:10.1103/PhysRevB.59.1758. [194] J.P. Perdew, J.A. Chevary, S.H. Vosko, K.A. Jackson, M.R. Pederson, D.J. Singh, C. Fiolhais, Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation, Phys. Rev. B. 46 (1992) 6671–6687. doi:10.1103/PhysRevB.46.6671. [195] J.P.R. De Villiers, Crystal structures of aragonite , strontianite , and witherite, Am. Mineral. 56 (1971) 758–767. http://www.minsocam.org/ammin/AM56/AM56_758.pdf. [196] W.H. Press, Numerical Recipes 3rd Edition: The Art of Scientific Computing, Cambridge Univ. Press, 2007. doi:10.2307/1269484. [197] H. Hellmann, D. Andrae, Hans Hellmann: Einführung in die Quantenchemie: Mit biografischen Notizen von Hans Hellmann jr., Springer-Verlag, 2015. [198] J. H., A Climbing Image Nudged Elastic Band Method for Finding Saddle Points and 9901. Minimum http://scitation.aip.org/content/aip/journal/jcp/113/22/10.1063/1.1329672 (accessed May 25, 2016). Energy (2000) Phys. 113 Paths, J. Chem. 147 [199] V. Fiorentini, M. Methfessel, Extracting convergent surface energies from slab calculations, J. Phys. Condens. Matter. 8 (1996) 6525–6529. doi:10.1088/0953-8984/8/36/005. [200] L. Li, C. Ortiz, Pervasive nanoscale deformation twinning as a catalyst for efficient energy dissipation in a bioceramic armour., Nat. Mater. 13 (2014) 501–7. doi:10.1038/nmat3920. [201] J.A. Chisholm, P.D. Bristowe, A first principles investigation of stacking fault energies and bonding in wurtzite materials, J. Phys. Condens. Matter. 11 (1999) 5057–5063. doi:10.1088/0953-8984/11/26/308. [202] W. Sun, V. Ageh, H. Mohseni, T.W. Scharf, J. Du, Experimental and computational studies (2014) 241903. titanate, Appl. Phys. Lett. 104 in zinc on stacking doi:10.1063/1.4883747. faults [203] S. Yoo, K.T. Butler, A. Soon, A. Abbas, J.M. Walls, A. Walsh, Identification of critical stacking faults in thin-film CdTe solar cells, Appl. Phys. Lett. 105 (2014) 062104. doi:10.1063/1.4892844. [204] P. Lazar, R. Podloucky, Mechanical properties of superhard BC5, Appl. Phys. Lett. 94 (2009) 251904. doi:10.1063/1.3159627. [205] J. a. Chisholm, P.D. Bristowe, Stacking fault energies in Si doped GaN: A first principles study, Appl. Phys. Lett. 77 (2000) 534. doi:10.1063/1.127035. [206] E.B. Tadmor, N. Bernstein, A first-principles measure for the twinnability of FCC metals, J. Mech. Phys. Solids. 52 (2004) 2507–2519. doi:10.1016/j.jmps.2004.05.002. [207] J.R. Rice, Dislocation nucleation from a crack tip: An analysis based on the Peierls concept, J. Mech. Phys. Solids. 40 (1992) 239–271. doi:10.1016/S0022-5096(05)80012-2. [208] R.J. Asaro, S. Suresh, Mechanistic models for the activation volume and rate sensitivity in metals with nanocrystalline grains and nano-scale twins, Acta Mater. 53 (2005) 3369–3382. doi:10.1016/j.actamat.2005.03.047. [209] Y.T. Zhu, X.Z. Liao, X.L. Wu, Deformation twinning in nanocrystalline materials, Prog. Mater. Sci. 57 (2012) 1–62. doi:10.1016/j.pmatsci.2011.05.001. [210] W. Li, Doctoral Thesis: First-principles description of planar faults in metals and alloys, http://www.diva- KTH portal.org/smash/record.jsf?pid=diva2:759534 (accessed May 25, 2016). Technology, Institute Royal of 2014. [211] Y. Juan, E. Kaxiras, Generalized stacking fault energy surfaces and dislocation properties of silicon: A first-principles theoretical study, Philos. Mag. A. 74 (1996) 1367–1384. doi:10.1080/01418619608240729. [212] P. Hirel, P. Marton, M. Mrovec, C. Elsässer, Theoretical investigation of {110} generalized stacking faults and their relation to dislocation behavior in perovskite oxides, Acta Mater. 58 (2010) 6072–6079. doi:10.1016/j.actamat.2010.07.025. 148 [213] M. Muzyk, Z. Pakiela, K.J. Kurzydlowski, Generalized stacking fault energy in magnesium theory calculations, Scr. Mater. 66 (2012) 219–222. alloys: Density functional doi:10.1016/j.scriptamat.2011.10.038. [214] X. Wu, R. Wang, S. Wang, Generalized-stacking-fault energy and surface properties for HCP metals: A first-principles study, Appl. Surf. Sci. 256 (2010) 3409–3412. doi:10.1016/j.apsusc.2009.12.042. [215] M.Y. Chou, M.L. Cohen, S.G. Louie, Theoretical study of stacking faults in silicon, Phys. Rev. B. 32 (1985) 7979–7987. doi:10.1103/PhysRevB.32.7979. [216] A.G. Crocker, The crystallography of deformation twinning in alpha-uranium, J. Nucl. Mater. 16 (1965) 306–326. doi:http://dx.doi.org/10.1016/0022-3115(65)90119-4. [217] J. Campbell, Castings, 2nd ed., Elsevier Butterworth-Heinemann, 2003. [218] A.M. Lovatt, D. Bassetti, H.R. Shercliff, Y. Bréchet, Process and alloy selection for 211–225. J. Cast Met. Res. (2000) 12 aluminium doi:10.1080/13640461.2000.11819358. casting, Int. [219] A. Lee, Y. Lu, A. Roche, T.Y. Pan, Influence of nano-structured silanols on the microstructure and mechanical properties of A4047 and A359 aluminum casting alloys, Int. J. Met. 10 (2016) 322–328. doi:10.1007/s40962-016-0044-4. [220] G.K. Sigworth, The Modification of Al-Si Casting Alloys: Important Practical and Theoretical Aspects, Int. J. Met. 2 (2008) 19–40. doi:10.1007/BF03355425. [221] R. Gopalan, N.K. Prabhu, Oxide bifilms in aluminium alloy castings – a review, Mater. Sci. Technol. 27 (2011) 1757–1769. doi:10.1179/1743284711Y.0000000033. [222] X. Cao, J. Campbell, Oxide inclusion defects in Al-Si-Mg cast alloys, Can. Metall. Q. 44 (2005) 435–448. doi:10.1179/cmq.2005.44.4.435. [223] J. Campbell, Melting, Remelting, and Casting for Clean Steel, Steel Res. Int. 88 (2017) 1600093. doi:10.1002/srin.201600093. [224] J. Campbell, The consolidation of metals: the origin of bifilms, J. Mater. Sci. 51 (2016) 96– 106. doi:10.1007/s10853-015-9399-9. [225] J. Campbell, An overview of the effects of bifilms on the structure and properties of cast alloys, Metall. Mater. Trans. B Process Metall. Mater. Process. Sci. 37 (2006) 857–863. doi:10.1007/BF02735006. [226] L. Liu, F.H. Samuel, Effect of inclusions on the tensile properties of Al-7% Si-0.35% Mg (1998) 2269–2281. J. Mater. Sci. 33 (A356.2) aluminium casting alloy, doi:10.1023/A:1004331219406. [227] J. Campbell, “Stop Pouring, Start Casting,” Int. J. Met. 6 (2012) 7–18. 149 doi:10.1007/BF03355529. [228] C. Reilly, N.R. Green, M.R. Jolly, The present state of modeling entrainment defects in the 611–628. process, Appl. Math. Model. (2013) shape doi:10.1016/J.APM.2012.04.032. casting 37 [229] C. Nyahumwa, N.R. Green, J.Campbell, Effect of mold-filling turbulence on fatigue properties of cast aluminum alloys, Trans. Am. Foundrymen’s Soc. 106 (1998) 215–223. https://www.scopus.com/record/display.uri?eid=2-s2.0- 0000104745&origin=inward&txGid=8457a50c577b2fc735fadc7b1b79fb22. [230] K. Wefers, C. Misra, Oxides and Hydroxides of Aluminum, Alcoa Tech. Pap. 19 (1987) 1– 100. http://epsc511.wustl.edu/Aluminum_Oxides_Alcoa1987.pdf. [231] M.A. Trunov, M. Schoenitz, E.L. Dreizin, Effect of polymorphic phase transformations in alumina layer on ignition of aluminium particles, Combust. Theory Model. 10 (2006) 603– 623. doi:10.1080/13647830600578506. [232] S.K. Lee, S.B. Lee, S.Y. Park, Y.S. Yi, C.W. Ahn, Structure of Amorphous Aluminum Oxide, Phys. Rev. Lett. 103 (2009) 95501–95504. doi:10.1103/PhysRevLett.103.095501. [233] S. Davis, G. Gutierrez, Structural, elastic, vibrational and electronic properties of amorphous Al2O3 from ab initio calculations, J Phys Condens Matter. 23 (2011) 495401. doi:10.1088/0953-8984/23/49/495401. [234] T.C.C. Chou, D. Adamson, T.G.G. Nieh, J. Mardinly, Microstructural evolution and properties of nanocrystalline alumina made by reactive sputtering deposition, Thin Solid Films. 205 (1991) 131–139. doi:10.1016/0040-6090(91)90294-8. [235] L.P.H. Jeurgens, W.G. Sloof, F.D. Tichelaar, E.J. Mittemeijer, Growth kinetics and mechanisms of aluminum-oxide films formed by thermal oxidation of aluminum, J. Appl. Phys. 92 (2002) 1649–1656. doi:10.1063/1.1491591. [236] R. Raiszadeh, W.D. Griffiths, A method to study the history of a double oxide film defect in liquid aluminum alloys, Metall. Mater. Trans. B Process Metall. Mater. Process. Sci. 37 (2006) 865–871. doi:10.1007/BF02735007. [237] W.D. Griffiths, R. Raiszadeh, Hydrogen, porosity and oxide film defects in liquid Al, J. Mater. Sci. 44 (2009) 3402–3407. doi:10.1007/s10853-009-3450-7. [238] S. Fox, J. Campbell, Visualisation of oxide film defects during solidification of aluminium alloys, Scr. Mater. 43 (2000) 881–886. doi:10.1016/S1359-6462(00)00506-6. [239] Q.G. Wang, D. Apelian, D.A. Lados, Fatigue behavior of A356-T6 aluminum cast alloys. Part I. Effect of casting defects, J. Light Met. 1 (2001) 73–84. doi:10.1016/S1471- 5317(00)00008-0. [240] A. Mitrasinovic, F.C. Robles Hernández, M. Djurdjevic, J.H. Sokolowski, On-line 150 prediction of the melt hydrogen and casting porosity level in 319 aluminum alloy using thermal analysis, Mater. Sci. Eng. A. 428 (2006) 41–46. doi:10.1016/j.msea.2006.04.084. [241] E. Calvié, J. Réthoré, L. Joly-Pottuz, S. Meille, J. Chevalier, V. Garnier, Y. Jorand, C. Esnouf, T. Epicier, J.B. Quirk, K. Masenelli-Varlot, Mechanical behavior law of ceramic nanoparticles from transmission electron microscopy in situ nano-compression tests, Mater. Lett. 119 (2014) 107–110. doi:10.1016/j.matlet.2014.01.002. [242] E. Calvié, L. Joly-Pottuz, C. Esnouf, P. Clément, V. Garnier, J. Chevalier, Y. Jorand, A. Malchère, T. Epicier, K. Masenelli-Varlot, Real time TEM observation of alumina ceramic nano-particles during compression, J. Eur. Ceram. Soc. 32 (2012) 2067–2071. doi:10.1016/j.jeurceramsoc.2012.02.029. [243] K. Chenoweth, A.C.T. van Duin, W.A. Goddard, ReaxFF Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation, J. Phys. Chem. A. 112 (2008) 1040–1053. doi:10.1021/jp709896w. [244] M.P. Tosi, F.G. Fumi, Ionic sizes and born repulsive parameters in the NaCl-type alkali halides—II: The generalized Huggins-Mayer form, J. Phys. Chem. Solids. 25 (1964) 45– 52. doi:10.1016/0022-3697(64)90160-X. [245] K. Muralidharan, K.-D. Oh, P.A. Deymier, K. Runge, J.H. Simmons, Molecular dynamics simulations of atomic-level brittle fracture mechanisms in amorphous silica, J. Mater. Sci. 42 (2007) 4159–4169. doi:10.1007/s10853-007-1638-2. [246] B. Vessal, M. Amini, C.R.A. Catlow, Computer simulation of the structure of silica glass, J. Non. Cryst. Solids. 159 (1993) 184–186. doi:10.1016/0022-3093(93)91295-E. [247] B.P. Feuston, S.H. Garofalini, Empirical three‐body potential for vitreous silica, J. Chem. Phys. 89 (1988) 5818–5824. doi:10.1063/1.455531. [248] S.K. Mitra, M. Amini, D. Fincham, R.W. Hockney, Molecular dynamics simulation of 365–372. Philos. Mag. (1981) 43 silicon doi:10.1080/13642818108221906. dioxide glass, B. [249] T.P. Senftle, S. Hong, M.M. Islam, S.B. Kylasa, Y. Zheng, Y.K. Shin, C. Junkermeier, R. Engel-Herbert, M.J. Janik, H.M. Aktulga, T. Verstraelen, A. Grama, A.C.T. van Duin, The ReaxFF reactive force-field: development, applications and future directions, Npj Comput. Mater. 2 (2016) 15011. doi:10.1038/npjcompumats.2015.11. [250] W.G. Hoover, Computational Statistical Mechanics, Elsevier, 1996. doi:10.1016/0097- 8485(92)80060-D. [251] H.M. Aktulga, J.C. Fogarty, S.A. Pandit, A.Y. Grama, Parallel reactive molecular dynamics: Numerical methods and algorithmic techniques, Parallel Comput. 38 (2012) 245–259. doi:10.1016/j.parco.2011.08.005. [252] S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, J. Comput. 151 Phys. 117 (1995) 1–19. doi:10.1006/jcph.1995.1039. [253] J.R. Morris, C.Z. Wang, K.M. Ho, C.T. Chan, Melting line of aluminum from simulations of coexisting phases, Phys. Rev. B. 49 (1994) 3109–3115. doi:10.1103/PhysRevB.49.3109. [254] X. Shi-Gang, S. Li-Xin, Z. Rong-Gen, H. Xing-Fang, Properties of aluminium oxide coating on aluminium alloy produced by micro-arc oxidation, Surf. Coatings Technol. 199 (2005) 184–188. doi:10.1016/j.surfcoat.2004.11.044. [255] R. Kniep, P. Lamparter, S. Steeb, Structure of Anodic Oxide Coatings on Aluminum, Angew. Chemie. 101 (1989) 975–977. doi:10.1002/ange.19891010737. [256] C.J. ENGBERG, E.H. ZEHMS, Thermal Expansion of Al2O3, BeO, MgO, B4C, SiC, and TiC Above 1000° C., J. Am. Ceram. Soc. 42 (1959) 300–305. doi:10.1111/j.1151- 2916.1959.tb12958.x. [257] C. Landron, L. Hennet, T.E. Jenkins, G.N. Greaves, J.P. Coutures, A.K. Soper, Liquid alumina: Detailed atomic coordination determined from neutron diffraction data using empirical potential structure refinement, Phys. Rev. Lett. 86 (2001) 4839–4842. doi:10.1103/PhysRevLett.86.4839. [258] G.T. Mearini, R.W. Hoffman, Tensile properties of aluminum/alumina multi-layered thin films, J. Electron. Mater. 22 (1993) 623–629. doi:10.1007/BF02666408. [259] Q. Zhang, X. Xiao, Y.T. Cheng, M.W. Verbrugge, A non-destructive method for measuring the mechanical properties of ultrathin films prepared by atomic layer deposition, Appl. Phys. Lett. 105 (2014) 061901. doi:10.1063/1.4892539. [260] R.G. Munro, Evaluated material properties for a sintered alpha-alumina, J. Am. Ceram. Soc. 80 (1997) 1919–1928. doi:10.1111/j.1151-2916.1997.tb03074.x. [261] D.H. Chung, G. Simmons, Pressure and Temperature Dependences of the Isotropic Elastic (1968) 5316–5326. J. Appl. Phys. 39 Moduli of Polycrystalline Alumina, doi:10.1063/1.1655961. [262] S.L. Shang, Y. Wang, Z.K. Liu, First-principles elastic constants of alpha- and theta-Al2O3, Appl. Phys. Lett. 90 (2007) 101909. doi:10190910.1063/1.2711762. [263] M.K. Tripp, C. Stampfer, D.C. Miller, T. Helbling, C.F. Herrmann, C. Hierold, K. Gall, S.M. George, V.M. Bright, The mechanical properties of atomic layer deposited alumina for use in micro- and nano-electromechanical systems, Sensors Actuators, A Phys. 130–131 (2006) 419–429. doi:10.1016/j.sna.2006.01.029. [264] K. Tapily, J.E. Jakes, D.S. Stone, P. Shrestha, D. Gu, H. Baumgart, A.A. Elmustafa, Nanoindentation Investigation of HfO2 and Al2O3 Films Grown by Atomic Layer Deposition, J. Electrochem. Soc. 155 (2008) H545. doi:10.1149/1.2919106. [265] C.F. Herrmann, F.W. DelRio, S.M. George, V.M. Bright, Properties of atomic layer 152 deposited Al2O3/ZnO dielectric films grown at low temperature for RF MEMS, Micromach. Microfabr. Process Technol. X. 5715 (2005) 159–166. doi:10.1117/12.589322. [266] G. Alcalá, P. Skeldon, G.E. Thompson, a B. Mann, H. Habazaki, K. Shimizu, Mechanical properties of amorphous anodic alumina and tantala films using nanoindetation, Nanotechnology. 13 (2002) 451–455. doi:10.1088/0957-4484/13/4/302. [267] Z. Xia, L. Riester, B.W. Sheldon, W. a Curtin, J. Liang, A. Yin, J.M. Xu, Mechanical properties of highly ordered nanoporous anodic alumina membranes, Rev. Adv. Mater. Sci. 6 (2004) 131–139. http://www.ipme.ru/e-journals/RAMS/no_2604/sheldon/sheldon.pdf (accessed August 5, 2018). [268] M.R. Gallas, G.J. Piermarini, Bulk Modulus and Young’s Modulus of Nanocrystalline γ- (1994) 2917–2920. doi:10.1111/j.1151- J. Am. Ceram. Soc. 77 Alumina, 2916.1994.tb04524.x. [269] G. Gutiérrez, B. Johansson, Molecular dynamics study of structural properties of amorphous Al2O3, Phys. Rev. B. 65 (2002) 104202. doi:10.1103/PhysRevB.65.104202. [270] J.B. Wachtman, T.G. Scuderi, G.W. Cleek, Linear Thermal Expansion of Aluminum Oxide and Thorium Oxide from 100° to 1100° K, J. Am. Ceram. Soc. 45 (1962) 319–323. doi:10.1111/j.1151-2916.1962.tb11159.x. [271] G. Kresse, M. Schmid, E. Napetschnig, M. Shishkin, L. Köhler, P. Varga, Structure of the Ultrathin Aluminum Oxide Film on NiAl(110), Science. 308 (2005) 1440–1442. doi:10.1126/science.1107783. [272] J. Mi, R.A. Harding, J. Campbell, Effects of the entrained surface film on the reliability of castings, Metall. Mater. Trans. A. 35 (2004) 2893–2902. doi:10.1007/s11661-004-0237-y. [273] X. Cao, J. Campbell, The nucleation of Fe-Rich phases on oxide films in Al-11.5Si-0.4Mg cast alloys, Metall. Mater. Trans. A. 34 (2003) 1409–1420. doi:10.1007/s11661-003-0253- 3. [274] J. Karch, R. Birringer, H. Gleiter, Ceramics ductile at low temperature, Nature. 330 (1987) 556–558. doi:10.1038/330556a0. [275] D. Dispinar, Determination of metal quality of aluminium and its alloys, University of Birmingham, 2006. http://etheses.bham.ac.uk/10/ (accessed June 26, 2018). [276] J. Campbell, Entrainment defects, Mater. Sci. Technol. 22 (2006) 127–145. doi:10.1179/174328406X74248. [277] A. Needleman, A Continuum Model for Void Nucleation by Inclusion Debonding, J. Appl. Mech. 54 (1987) 525. doi:10.1115/1.3173064. [278] J. Liu, Z. Huang, Z. Pan, Q. Wei, X. Li, Y. Qi, Atomistic Origin of Deformation Twinning 105501. Biomineral Aragonite, (2017) Lett. 118 in Phys. Rev. 153 doi:10.1103/PhysRevLett.118.105501. [279] S. Xia, Y. Qi, T. Perry, K.S. Kim, Strength characterization of Al/Si interfaces: A hybrid method of nanoindentation and finite element analysis, Acta Mater. 57 (2009) 695–707. doi:10.1016/j.actamat.2008.10.011. [280] M. Divandari, J. Campbell, Oxide film characteristics of Al–7Si–Mg alloy in dynamic 182–187. J. Cast Met. Res. (2004) 17 conditions doi:10.1179/136404604225017546. casting, Int. in [281] C. Nyahumwa, Multiple defect distributions on weibull statistical analysis of fatigue life of cast aluminium alloys, African J. Sci. Technol. 6 (2010) 43–54. doi:10.4314/ajst.v6i2.55174. [282] Q.G. Wang, P.E. Jones, Prediction of Fatigue Performance in Aluminum Shape Castings Containing Defects, Metall. Mater. Trans. B. 38 (2007) 615–621. doi:10.1007/s11663-007- 9051-4. [283] M. Armand, J.-M.J.-M. Tarascon, Building better batteries, Nature. 451 (2008) 652–657. doi:10.1038/451652a. [284] S. Chu, A. Majumdar, Opportunities and challenges for a sustainable energy future, Nature. 488 (2012) 294–303. doi:10.1038/nature11475. [285] J. Janek, W.G. Zeier, A solid future for battery development, Nat. Energy. 1 (2016) 16141. doi:10.1038/nenergy.2016.141. [286] X. Han, Y. Gong, K. (Kelvin) Fu, X. He, G.T. Hitz, J. Dai, A. Pearse, B. Liu, H. Wang, G. Rubloff, Y. Mo, V. Thangadurai, E.D. Wachsman, L. Hu, Negating interfacial impedance in garnet-based solid-state Li metal batteries, Nat. Mater. 16 (2016) 572–579. doi:10.1038/nmat4821. [287] W. Luo, Y. Gong, Y. Zhu, K.K. Fu, J. Dai, S.D. Lacey, C. Wang, B. Liu, X. Han, Y. Mo, E.D. Wachsman, L. Hu, Transition from Superlithiophobicity to Superlithiophilicity of Garnet Solid-State Electrolyte, J. Am. Chem. Soc. 138 (2016) 12258–12262. doi:10.1021/jacs.6b06777. [288] K. (Kelvin) Fu, Y. Gong, B. Liu, Y. Zhu, S. Xu, Y. Yao, W. Luo, C. Wang, S.D. Lacey, J. Dai, Y. Chen, Y. Mo, E. Wachsman, L. Hu, Toward garnet electrolyte–based Li metal batteries: An ultrathin, highly effective, artificial solid-state electrolyte/metallic Li interface, Sci. Adv. 3 (2017) e1601659. doi:10.1126/sciadv.1601659. [289] D. Lin, Y. Liu, Z. Liang, H.-W. Lee, J. Sun, H. Wang, K. Yan, J. Xie, Y. Cui, Layered reduced graphene oxide with nanoscale interlayer gaps as a stable host for lithium metal anodes, Nat. Nanotechnol. 11 (2016) 626–632. doi:10.1038/nnano.2016.32. [290] G.M. Stone, S.A. Mullin, A.A. Teran, D.T. Hallinan, A.M. Minor, A. Hexemer, N.P. Balsara, Resolution of the Modulus versus Adhesion Dilemma in Solid Polymer Electrolytes for Rechargeable Lithium Metal Batteries, J. Electrochem. Soc. 159 (2012) 154 A222–A227. doi:10.1149/2.030203jes. [291] Y. TAKEDA, O. YAMAMOTO, N. IMANISHI, Lithium Dendrite Formation on a Lithium Metal Anode from Liquid, Polymer and Solid Electrolytes, Electrochemistry. 84 (2016) 210–218. doi:10.5796/electrochemistry.84.210. [292] J. Wandt, C. Marino, H.A. Gasteiger, P. Jakes, R.-A. Eichel, J. Granwehr, Operando electron paramagnetic resonance spectroscopy – formation of mossy lithium on lithium anodes during charge–discharge cycling, Energy Environ. Sci. 8 (2015) 1358–1367. doi:10.1039/C4EE02730B. [293] J. Fang, H. You, C. Zhu, P. Kong, M. Shi, X. Song, B. Ding, Thermodynamic and kinetic competition in silver dendrite growth, Chem. Phys. Lett. 439 (2007) 204–208. doi:10.1016/J.CPLETT.2007.03.046. [294] D. Aurbach, E. Zinigrad, Y. Cohen, H. Teller, A short review of failure mechanisms of lithium metal and lithiated graphite anodes in liquid electrolyte solutions, Solid State Ionics. 148 (2002) 405–416. doi:10.1016/S0167-2738(02)00080-2. [295] F. Yonemoto, A. Nishimura, M. Motoyama, N. Tsuchimine, S. Kobayashi, Y. Iriyama, Temperature effects on cycling stability of Li plating/stripping on Ta-doped Li7La3Zr2O12, J. Power Sources. 343 (2017) 207–215. doi:10.1016/J.JPOWSOUR.2017.01.009. [296] N. Nitta, F. Wu, J.T. Lee, G. Yushin, Li-ion battery materials: present and future, Mater. Today. 18 (2015) 252–264. doi:10.1016/J.MATTOD.2014.10.040. [297] M.D. Radin, J.F. Rodriguez, F. Tian, D.J. Siegel, Lithium Peroxide Surfaces Are Metallic, While Lithium Oxide Surfaces Are Not, J. Am. Chem. Soc. 134 (2012) 1093–1103. doi:10.1021/ja208944x. [298] J.S. Gnanaraj, R.W. Thompson, S.N. Iaconatti, J.F. DiCarlo, K.M. Abraham, Formation and Growth of Surface Films on Graphitic Anode Materials for Li-Ion Batteries, Electrochem. Solid-State Lett. 8 (2005) A128. doi:10.1149/1.1850390. [299] S.-K. Jeong, M. Inaba, R. Mogi, Y. Iriyama, T. Abe, Z. Ogumi, Surface Film Formation on a Graphite Negative Electrode in Lithium-Ion Batteries: Atomic Force Microscopy Study on the Effects of Film-Forming Additives in Propylene Carbonate Solutions, Langmuir. 17 (2001) 8281–8286. doi:10.1021/la015553h. [300] V.A. Agubra, J.W. Fergus, The formation and stability of the solid electrolyte interface on 153–162. Sources. (2014) Power 268 graphite the doi:10.1016/J.JPOWSOUR.2014.06.024. anode, J. [301] A.A. Talin, D. Ruzmetov, A. Kolmakov, K. McKelvey, N. Ware, F. El Gabaly, B. Dunn, H.S. White, Fabrication, Testing, and Simulation of All-Solid-State Three-Dimensional Li- Ion Batteries, ACS Appl. Mater. 32385–32391. doi:10.1021/acsami.6b12244. Interfaces. (2016) 8 155 [302] L. Oniciu, L. Mureşan, Some fundamental aspects of levelling and brightening in metal electrodeposition, J. Appl. Electrochem. 21 (1991) 565–574. doi:10.1007/BF01024843. [303] A. Hasnaoui, O. Politano, J.M.M. Salazar, G. Aral, R.K.K. Kalia, A. Nakano, P. Vashishta, Molecular dynamics simulations of the nano-scale room-temperature oxidation of aluminum single crystals, Surf. Sci. 579 (2005) 47–57. doi:10.1016/j.susc.2005.01.043. [304] B.E. Deal, A.S. Grove, General relationship for the thermal oxidation of silicon, J. Appl. Phys. 36 (1965) 3770–3778. doi:10.1063/1.1713945. [305] E. Dologlou, Self-diffusion in solid lithium, Glas. Phys. Chem. 36 (2010) 570–574. doi:10.1134/S1087659610050056. [306] Y. Oishi, Y. Kamei, M. Akiyama, T. Yanagi, Self-diffusion coefficient of lithium in lithium oxide, J. Nucl. Mater. 87 (1979) 341–344. doi:10.1016/0022-3115(79)90570-1. [307] H. Zheng, Y. Liu, S.X. Mao, J. Wang, J.Y. Huang, Beam-assisted large elongation of in situ formed Li2O nanowires, Sci. Reports 2012 2. 2 (2012) 542. doi:10.1038/srep00542. [308] P. Jagtap, A. Chakraborty, P. Eisenlohr, P. Kumar, Identification of whisker grain in Sn coatings by analyzing crystallographic micro-texture using electron back-scatter diffraction, Acta Mater. 134 (2017) 346–359. doi:10.1016/J.ACTAMAT.2017.05.063. [309] J. Belak, J.N. Glosli, D.B. Boercker, I.F. Stowers, Molecular Dynamics Simulation of Mechanical Deformation of Ultra-Thin Metal and Ceramic Films, MRS Proc. 389 (1995) 181. doi:10.1557/PROC-389-181. [310] I. Salehinia, J. Wang, D.F. Bahr, H.M. Zbib, Molecular dynamics simulations of plastic (2014) 119–132. in Nb/NbC multilayers, J. Plast. 59 Int. deformation doi:10.1016/J.IJPLAS.2014.03.010. [311] D. Wolf, V. Yamakov, S.R. Phillpot, A. Mukherjee, H. Gleiter, Deformation of nanocrystalline materials by molecular-dynamics simulation: relationship to experiments?, Acta Mater. 53 (2005) 1–40. doi:10.1016/J.ACTAMAT.2004.08.045. [312] P.S. Branício, J.-P. Rino, Large deformation and amorphization of Ni nanowires under uniaxial strain: A molecular dynamics study, Phys. Rev. B. 62 (2000) 16950–16955. doi:10.1103/PhysRevB.62.16950. [313] M.M. Islam, A. Ostadhossein, O. Borodin, A.T. Yeates, W.W. Tipton, R.G. Hennig, N. Kumar, A.C.T. van Duin, ReaxFF molecular dynamics simulations on lithiated sulfur cathode materials, 3383–3393. doi:10.1039/C4CP04532G. Phys. Chem. Chem. (2015) Phys. 17 [314] L. Van Brutzel, C.L. Rountree, R.K. Kalia, A. Nakano, P. Vashishta, Dynamic Fracture Mechanisms in Nanostructured and Amorphous Silica Glasses Million-Atom Molecular Dynamics Simulations, MRS Proc. 703 (2001) V3.9. doi:10.1557/PROC-703-V3.9. 156 [315] F. Yuan, L. Huang, Brittle to Ductile Transition in Densified Silica Glass, Sci. Rep. 4 (2015) 5035. doi:10.1038/srep05035. [316] R.T. Sanderson, An Interpretation of Bond Lengths and a Classification of Bonds, Science (80-. ). 114 (1951) 670–672. doi:10.1126/science.114.2973.670. [317] W.J. Mortier, S.K. Ghosh, S. Shankar, Electronegativity-equalization method for the calculation of atomic charges in molecules, J. Am. Chem. Soc. 108 (1986) 4315–4320. doi:10.1021/ja00275a013. [318] W.J. Mortier, K. Van Genechten, J. Gasteiger, Electronegativity equalization: application and parametrization, J. Am. Chem. Soc. 107 (1985) 829–835. doi:10.1021/ja00290a017. [319] A.K. Rappé, W.A. Goddard III, Charge Equilibration for Molecular Dynamics Simulations, J. Phys. Chem. 95 (1991) 3358–3363. doi:10.1021/j100161a070. [320] B. Mortazavi, A. Ostadhossein, T. Rabczuk, A.C.T. van Duin, Mechanical response of all- MoS 2 single-layer heterostructures: a ReaxFF investigation, Phys. Chem. Chem. Phys. 18 (2016) 23695–23701. doi:10.1039/C6CP03612K. [321] G.M. Odegard, B.D. Jensen, S. Gowtham, J. Wu, J. He, Z. Zhang, Predicting mechanical response of crosslinked epoxy using ReaxFF, Chem. Phys. Lett. 591 (2014) 175–178. doi:10.1016/J.CPLETT.2013.11.036. [322] B.D. Jensen, K.E. Wise, G.M. Odegard, The effect of time step, thermostat, and strain rate on ReaxFF simulations of mechanical failure in diamond, graphene, and carbon nanotube, J. Comput. Chem. 36 (2015) 1587–1596. doi:10.1002/jcc.23970. [323] A. Ostadhossein, S.-Y. Kim, E.D. Cubuk, Y. Qi, A.C.T. van Duin, Atomic Insight into the Lithium Storage and Diffusion Mechanism of SiO 2 /Al 2 O 3 Electrodes of Lithium Ion Batteries: ReaxFF Reactive Force Field Modeling, J. Phys. Chem. A. 120 (2016) 2114– 2127. doi:10.1021/acs.jpca.5b11908. [324] M. Biology, T. Scripps, CHARMM Fluctuating Charge Force Field for Proteins : I Parameterization and Application to Bulk Organic, J. Comput. Chem. 25 (2004) 1–16. http://onlinelibrary.wiley.com/doi/10.1002/jcc.10355/full (accessed November 13, 2016). [325] R. Chelli, P. Procacci, R. Righini, S. Califano, Electrical response in chemical potential equalization schemes, J. Chem. Phys. 111 (1999) 8569. doi:10.1063/1.480198. [326] R.A. Nistor, J.G. Polihronov, M.H. Müser, N.J. Mosey, A generalization of the charge (2006). for nonmetallic materials, J. Chem. Phys. 125 equilibration method doi:10.1063/1.2346671. [327] W.B. Dapp, M.H. Müser, Redox reactions with empirical potentials: Atomistic battery discharge simulations, J. Chem. Phys. 139 (2013). doi:10.1063/1.4817772. [328] N. Masaki, K. Doi, S. Nasu, T. Tanifuji, K. Uchida, Gitterstruktur der Oxide, Sulfide, 157 Selenide und Telluride des Lithiums, Natriums und Kaliums, J. Nucl. Mater. 84 (1979) 341– 342. [329] C. GANDHI, M.F. ASHBY, FRACTURE-MECHANISM MAPS FOR MATERIALS WHICH CLEAVE: F.C.C., B.C.C. AND H.C.P. METALS AND CERAMICS, in: Perspect. Creep Fract., Pergamon, 1983: pp. 33–70. doi:10.1016/b978-0-08-030541-7.50006-0. [330] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, Y. Bengio, Generative Adversarial Nets, in: Z. Ghahramani, M. Welling, C. Cortes, N.D. Lawrence, K.Q. Weinberger (Eds.), Adv. Neural Inf. Process. Syst. 27, Curran Associates, Inc., 2014: pp. 2672–2680. http://papers.nips.cc/paper/5423-generative-adversarial- nets.pdf. [331] P. Isola, J.-Y. Zhu, T. Zhou, A.A. Efros, Image-To-Image Translation With Conditional Adversarial Networks, in: IEEE Conf. Comput. Vis. Pattern Recognit., 2017. [332] B.D. Tracey, K. Duraisamy, J.J. Alonso, A Machine Learning Strategy to Assist Turbulence Model Development, in: 53rd AIAA Aerosp. Sci. Meet., American Institute of Aeronautics and Astronautics, Reston, Virginia, 2015. doi:10.2514/6.2015-1287. [333] F. Brockherde, L. Vogt, L. Li, M.E. Tuckerman, K. Burke, K.R. Müller, Bypassing the (2017). learning, Nat. Commun. 8 Kohn-Sham doi:10.1038/s41467-017-00839-3. equations with machine 158