Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials
(USC Thesis Other)
Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
QUANTUM MOLECULAR DYNAMICS AND MACHINE LEARNING FOR NON-EQUILIBRIUM PROCESSES IN QUANTUM MATERIALS By Subodh C. Tiwari _____________________________________________________________________ A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (MATERIALS SCIENCE) May 2020 Copyright 2020 Subodh C. Tiwari ii Acknowledgements First, I would like to express my sincere gratitude towards my Ph.D. advisor, Prof. Priya Vashishta for his continuous guidance and encouragement during last 6 years. I am extremely fortunate to have him as my advisor. His invaluable training has guided me on all aspect of my Ph.D. and helped me become a better researcher. I am very thankful to Prof. Aiichiro Nakano for his immense support both in my research and in my life. His hard work, dedication and discipline in research has always encouraged me to do better than yesterday. His modesty and caring personality always made me humble in my personal and professional life. I admired Prof. Rajiv Kalia for his love for fundamental physics and thought-provoking questions. I would like to also thank Prof. Paulo Branicio, Dr. Ken-ichi Nomura and Prof. Fuyuki Shimojo. Without their constant support in research and innovative way to approach most of my technical problem, it would not be possible to complete my research. I would also like to thank my colleagues in our research group (CACS), USC and Fukuoka University, who have made my journey during this whole endeavor extremely enjoyable. In CACS, there was always someone to hear and discuss our problems in research, whether it’s debate in the office or discussion over drinks. My friends at USC who always celebrated my success with me both big and small and supported me though my setbacks. Finally, I am extremely thankful to my family and friends in India for their unconditional love and support, and their constant encouragement in my scientific endeavors. I am especially grateful to my sister and my mother whose hard work allowed me to pursue my scientific aspirations. iii Table of Contents Acknowledgements ......................................................................................................................... ii List of Figures ................................................................................................................................ vi Abstract .......................................................................................................................................... xi 1. Introduction ............................................................................................................................. 1 1.1. Light driven non-thermal amorphization process in GeTe and Sb2Te3 .......................... 2 1.2. Hydrogen bond is key to shock resilience of Aramid fibers ........................................... 4 1.3. Multiple reaction pathways in shocked 2,4,6-triamino- 1,3,5 trinitrobenzene crystal.... 5 1.4. Electric field induced polarization switching in AlxSc(1-x)N ........................................... 6 1.5. Machine learning solution of Schrödinger equation via principle of least action .......... 7 2. Method .................................................................................................................................... 9 2.1. Overview of molecular dynamics ................................................................................... 9 2.2. Interatomic forcefield.................................................................................................... 10 2.3. Reactive forcefield ........................................................................................................ 11 2.4. Quantum molecular dynamics ...................................................................................... 13 2.5. Non-adiabatic quantum molecular dynamics................................................................ 16 2.5.1. Long-range corrected ground state dynamics ........................................................... 17 2.5.2. Electronic excitations ................................................................................................ 17 2.5.3. Force in molecular dynamics .................................................................................... 19 2.6. Multiscale shock Theory (MSST) ................................................................................. 22 2.7. Finite electric field in quantum molecular dynamics .................................................... 24 3. Light Driven Non-thermal Amorphization in Sb2Te3 and GeTe........................................... 26 3.1. Introduction ................................................................................................................... 26 3.2. Non-thermal amorphization of Sb2Te3.......................................................................... 28 3.2.1. Simulation setup........................................................................................................ 28 3.2.2. Analysis of crystalline to amorphous phase change ................................................. 28 3.2.3. Evolution of excited state.......................................................................................... 32 3.3. Non-thermal amrophization of GeTe ............................................................................ 36 3.3.1. Simulation setup........................................................................................................ 36 3.3.2. Amorphous to crystalline transformation ................................................................. 37 3.3.3. Electronic and ionic transformation during excited state dynamics ......................... 40 3.4. Conclusion .................................................................................................................... 43 4. Shock Simulation of Aramid Fibers ..................................................................................... 44 4.1. Introduction ................................................................................................................... 44 iv 4.2. MSST simulation with quantum molecular dynamics .................................................. 46 4.2.1. Simulation details...................................................................................................... 46 4.2.2. Hugoniot curve.......................................................................................................... 48 4.2.3. Structural transformation .......................................................................................... 49 4.2.4. Global structure search using active learning ........................................................... 53 4.2.5. Analysis of predicted structure ................................................................................. 56 4.3. Reactive forcefield modeling of Kevlar ........................................................................ 58 4.3.1. Development of ReaxFF-lg Forcefield ..................................................................... 59 4.3.2. Bond dissociations energy ........................................................................................ 60 4.3.3. Energy volume curve ................................................................................................ 63 4.3.4. Characterization of PPTA crystal ............................................................................. 64 4.4. MSST simulation with Reactive molecule dynamics ................................................... 66 4.4.1. Simulation setup........................................................................................................ 66 4.4.2. Finite size effect of simulation cell ........................................................................... 67 4.4.3. Shock Hugoniot ........................................................................................................ 68 4.4.4. Paracrystalline transformation under mild shock ..................................................... 69 4.5. Conclusion .................................................................................................................... 70 5. Multiple Reaction Pathway in Shocked 2,4,6-Triamino-1,3,5-Trinitrobenzene crystal ....... 71 5.1. Introduction ................................................................................................................... 71 5.2. Simulation details.......................................................................................................... 73 5.3. Results ........................................................................................................................... 74 5.3.1. Temperature profile .................................................................................................. 74 5.3.2. Fragment Analysis .................................................................................................... 75 5.3.3. Detection of new reaction pathway........................................................................... 80 5.4. Conclusion .................................................................................................................... 84 6. Electric field induced spontaneous polarization switching in Al(1-x)ScxN ............................ 85 6.1. Introduction ................................................................................................................... 85 6.2. Computational synthesis of Al(1-x)ScxN alloys .............................................................. 86 6.2.1. Local structure of Al(1-x)ScxN alloys ......................................................................... 86 6.2.2. Phase diagram of Al(1-x)ScxN alloys .......................................................................... 87 6.3. Switching pathway of Al(1-x)ScxN alloys ...................................................................... 88 6.3.1. Energy barrier calculation ......................................................................................... 88 6.3.2. Switching Metastable Layered nitride ...................................................................... 89 6.4. Conclusion .................................................................................................................... 91 7. Machine Learning Assisted Principle of Least Action Formalism for Schrodinger Equation..92 7.1. Introduction ................................................................................................................... 92 7.2. Method .......................................................................................................................... 94 7.3. Results ........................................................................................................................... 98 7.4. Conclusion .................................................................................................................. 101 v 8. Conclusion .......................................................................................................................... 102 References ................................................................................................................................... 104 vi List of Figures Figure 1: Scheme Set and Reset process in phase change materials a) via melt-quench process using laser b) via non-thermal photoexcitation. 3 Figure 2: hysteresis curves of Al1-xScxN (x = 0.27 – 0.43) 7 Figure 3: Dimensionless Lennard-Jones potential for inert gases. 11 Figure 4: Energy terms included in ReaxFF where Ebond is the bond energy, Elp is the lone pair energy, Eover and Eunder are the over- and under- coordination energies, Eval is the valence angle energy, Epen is the penalty energy, Ecoa is the three-body conjugation energy, Etors is the torsion angle energy, Econj is the four-body conjugation energy, Ehbond is the hydrogen bond interaction energy, EvdWaals, ECoulomb and Elg are the van der Waals, Coulomb energies and dispersion energy, respectively. 12 Figure 5: Workflow schematic for general QMD simulation to obtain Kohn-Sham eigenstates. After input of initial atomic positions {Ri} and atomic velocities {vi}, steps include making an initial guess for the charge density 𝜌 (r); computation of the system Hamiltonian, including kinetic energy T, Hartree potential VH, exchange-correlation potential VXC, local potential VL, and non-local potential VNL; solving the Kohn-Sham equations; computing new charge density 𝜌 ′ (r) checking convergence criterion; outputting new atomic positions if criterion is met or mixing ρ′(r) and 𝜌 (r) to begin to next iteration if not. 15 Figure 6: Sb2Te3 simulation supercell structures (a) initial crystalline structure of Sb2Te3, (b) structure at 300K (c-f) structure of Sb2Te3 after 4.22 ps excited state dynamics with n = 2.6, 5.2, 10.3 and 15.0 % excited valence electrons, respectively. Yellow and purple spheres correspond to antimony and tellurium atoms, respectively 29 Figure 7: Mean square displacement of Sb as a function of time. Cyan, green, orange and red curve corresponds to excitation of n=2.6%, n=5.2%, n= 10.3% and n=15% respectively at temperature 500 K. Black curve corresponds to the adiabatic MD at 500 K. n is percentage of excited electron out of the total electrons. 30 Figure 8: Analysis of the structural evolution (a) total radial distribution function as a function of number of excitons (b) evolution of the number of wrong bonds for different excitons (c) evolution of simulated diffraction pattern for n= 10.3% excitation. 31 Figure 9: Average Mulliken charge per component as a function of time. (a) n=2.6%, (b) n = 5.2%, (c) n = 10.3% excitation (c) n = 15% excitation. 33 Figure 10: Evolution of Mulliken bond overlap analysis as a function of time for a) Sb- Sb, b) Sb-Te, and c) Te-Te bond pairs for n = 2.6%. Similarly, d-f) correspond to Sb-Sb, Sb-Te and Te-Te bond pairs for n = 5.2%, g-i) correspond to Sb-Sb, Sb-Te and Te-Te bond pairs for n = 10.3%, and j-l) correspond to Sb-Sb, Sb-Te and Te-Te bond pairs for n = 15%, respectively. Black color curve corresponds to the t = -0.5 ps (before excitation). Red, green, cyan, and blue colors correspond to times t = 0 ps , 0.5 ps, 1.0 ps, and 2.5 ps after excitation. 35 vii Figure 11: GeTe simulation supercell structures (a) initial crystalline structure of Sb2Te3, (b) structure at 300K (c-f) structure of GeTe after 4.22 ps excited state dynamics with n = 2.6, 5.2, 10.3 and 15.0 % excited valence electrons, respectively. Orange and blue sphere represent germanium (Ge) and tellurium (Te) atoms, respectively 37 Figure 12: Mean square displacement of Ge atom type as a function of time. Cyan, green, orange and red curve corresponds to excitation of n=3.3%, n=5%, n= 7.5% and n=10% respectively at temperature 800 K. Black curve corresponds to the adiabatic MD at 800 K. n is percentage of excited electron out of the total electrons. 38 Figure 13: Analysis of the structural evolution (a) total radial distribution function as a function of number of excitons (b) evolution of the number of wrong bonds for different excitons (c) evolution of simulated diffraction pattern for n= 10.3% excitation. 39 Figure 14: Average mulliken charge per component as a function of time. (a) n=3.3%, (b) n = 10% valence band excitation excitation. 41 Figure 15: Evolution of Mulliken bond overlap analysis as a function of time for a) Ge- Ge, b) Ge-Te, and c) Te-Te bond pairs for n = 5%. Cyan color curve corresponds to the t = -0.1 ps (before excitation). orange, green curves correspond to times t = 0.1 ps , 2.0 ps after excitation. 5d) average mulliken bond overlap for Ge-Ge (red curve), Te-Te(blue curve) and Ge-Te(black curve) as a function of time. 42 Figure 16: Structure of aramid fibers and PPTA (p-phenylene terephthalamide) crystal. (a) Aramid fibers microstructure is composed of elongated, highly crystalline fibrils where PPTA polymer chains are aligned along the fibril length forming stacked PPTA sheets. (b) PPTA polymer chains are held together by in-plane hydrogen bonding (orange colored dotted lines) to create a PPTA sheet. Axes a, b, and c correspond to [100], [010], and [001] crystallographic directions. (c) PPTA crystal oriented along b, highlighting the PPTA sheet stacking structure. (d) PPTA crystal oriented along the polymer chain (fibril length). Grey, white, red, and blue spheres correspond to carbon, hydrogen, oxygen, and nitrogen atoms, respectively. 45 Figure 17: Shock Hugoniot along [100] (a) and [010] directions (b), respectively. PPTA structure in elastic (c), transformation (d), and cross-linking (e) regimes for shocks along [100] direction and figure (f-h) along [010] direction, respectively. Chains undergoing structural or trans-cis transformation are highlighted. Figure (i-k) Analysis of structural changes in different regimes by the evolution of dihedral angle between phenylene and amide groups in the elastic regime, fraction of hydrogen bonds preserved in the plastic/transformation regimes, and fraction of sp2 preserved in the cross-linking regime. 47 Figure 18: Shock-induced pressure in PPTA. (a) Induced pressure during shock along the [100] and (b) [010] directions. Red, blue, and green symbols indicate elastic, inelastic (transformation/plastic), and cross-linking regimes. 48 Figure 19: Dihedral angle between phenylene group and amide. We computed dihedral angle between the phenylene and amide groups as a function of time to observe the conformational changes during elastic shock 49 Figure 20: trans (a) and cis (b) conformations of a PPTA polymer chain. PPTA crystals have 100% trans conformation. 20% of polymer chains transform to cis during shock 51 viii along the [010] direction, leading to a planar amorphous phase. (c) Radial distribution function of polymer chains center of mass in the (001) plane for crystalline (C), planar amorphous (A), compressed crystalline (c-C) and transformed (T) systems. The insets show the corresponding center of mass configurations. Figure 21: (a) Energy landscape of the compressed system as a function of polymer chains’ reduced displacement in [010] and [001] directions. Reduced displacements are calculated using Gödel encoding (see Supporting information). (b) Predicted new phase of PPTA under shock along the [100] direction. All the phenylene groups in polymer chains are coplanar in this PPTA phase. (c) Unshocked PPTA monoclinic primitive unit cell. (d) Predicted new phase PPTA monoclinic primitive unit cell. 53 Figure 22: Bayesian optimization to find global minimum. 55 Figure 23: Energy difference of transformed state: We computed energy difference between transformed PPTA state and compressed PPTA state at different uniaxial compression in [100] direction. Transformed state shows lower energy for all used functionals. 56 Figure 24: Results of the ReaxFF refinement for the Kevlar systems: (a) the central C-N bond in the amide linkage before (right) and after (left) refining; (b) the C-N bond between the amide linkage and aromatic ring in the single Kevlar chain before (right) and after (left) refining; (c) the central H-N bond in the single Kevlar chain before (left) and after (right) refining. Note that energy contributions of the long-range-correction terms (Elg) were not included in the full bond dissociation curves above. 60 Figure 25: Bond dissociation energies described by the refined ReaxFF force field, DFT and ReaxFF-lg: (a) C-N single bond; (b) C-N double bond; (c) H-N single bond. Note that after the refinement, the ReaxFF force field still reproduces the C-N single/double bonds and H-N single bond, which were included in the previous training set. 61 Figure 26: Energy volume curve of Kevlar after planar compression in [100] and [010] 62 Figure 27: RMD simulations of thermal equilibrations for the Kevlar crystal structure: (a) at 10 K for 20 ps; (b) at 300 K for 125 ps. During the thermal process, the Kevlar structure kept the AB arrangements stable in the xy plane. In addition, no significant distortion or bending was observed in the amide linkage (the xz plane) for up to 300 K. 63 Figure 28: shock speed in the Kevlar after immediately Hugoniot shock limit as function of increasing size. 66 Figure 29: Shock Hugoniot along [100] (a) and [010] directions (b), respectively for RMD simulation. 67 Figure 30: Paracrystalline transformation of Kevlar during shock loading in [010] direction. 69 Figure 31: (a) Single TATB molecule where cyan, grey, red and blue represent carbon, hydrogen, oxygen and nitrogen, respectively. (b) TATB crystalline unit cell. Each unit cell consists of 2 TATB molecules. The unit cell has the lattice parameters of a = 9.010 Å, b = 9.028 Å, c = 6.812 Å, a = 108.58˚, b = 91.82˚ and g = 119.97˚.158 70 ix Figure 32: Schematic diagram of system setup. The orange-colored block at the left end of the system represents a repulsive flyer, while the right end is held by a reflective wall. 72 Figure 33: Temperature profile of a shocked TATB crystal at different times. A discontinuous temperature jump is observed between the unshocked (blue) and shocked (red) regions. 73 Figure 34: Temperature profile of the system as a function of the position along the a- axis and time. 74 Figure 35: Number of molecular fragments formed as a function of time for shock speed 9.98 km/s. Red, blue, black and green curves correspond to the formation of N2, H2O, NH3 and CO2, respectively. 75 Figure 36: (a) Fraction of N2 molecules that are formed by two distinct pathways as a function of time. The red line shows the fraction of N2 formed by N atoms from a single TATB molecule (intramolecular pathway). The blue line shows N2 originating from more than one TATB molecule (intermolecular pathway). (b) Fraction of H2O molecules that are formed by two distinct pathways as a function of time. The red line shows the fraction of H2O formed by H and O atoms from a single TATB molecule (intramolecular pathway). The blue line shows H2O originating from more than one TATB molecule (intermolecular pathway) 77 Figure 37: Intermediate species formed as a function of time. NO2, NO, HONO, OH and NH2 are represented by blue, green, red, black and magenta, respectively. 78 Figure 38: a 2D layer in the a-b crystallographic plane of TATB crystal. Cyan, grey, red and blue represent carbon, hydrogen, oxygen and nitrogen, respectively. Black circles enclose proximate intermolecular NO2-NH2 pairs 80 Figure 39: Time evolution of the fraction of N2 formed from different reaction pathways. 81 Figure 40: Time evolution of large clusters (the number of atoms in the cluster, n > 24). (a) Fraction of the number of atoms in large clusters to the total number of atoms. (b) Atomic composition of the large clusters as a function of time. Red and blue lines 83 Figure 41: local structure of Al1-xScxN ferroelectrics. a) Al0.67Sc0.33N local structure b) Al0.60Sc0.40N local structure. Cyan, black and brown spheres correspond to nitrogen, aluminum and scandium atoms. Red and blue sphere highlight ScN local structure. 85 Figure 42: Phase diagram of Al(1-x)ScxN alloys 86 Figure 43: Al0.75Sc0.25N configuration and energy barrier. Al0.75Sc0.25N in down (a)/up (b) position. (c) energy barrier for switching of Al0.75Sc0.25N (red curve) and Al0.75Sc0.25N (blue curve). Grey and yellow spheres correspond to Al and Sc atoms respectively. Blue/red sphere corresponds to Ndown/Nup. 88 Figure 44: position of Al (red curves) and Sc (blue curves) as a function of time for one layer. Snapshot of Sc and its surrounding N atoms at different time are shown by around image. Yellow and blue sphere corresponds to the Sc and N, respectively. 89 Figure 45: A scheme of feed forward neural network employed in our calculation. We have used tanh as a non-linear unit. 95 x Figure 46: Training results of neural network to obtain radial wave function using PLA method. a) loss as a function of epochs b-d) 1s, 2s, and 2p 97 Figure 47: Training results of neural network to obtain radial wave function using Hamiltonian method. a) loss as a function of epochs b-d) 1s, 2s, and 3s 98 Figure 48: SCF iteration and resulting 1s, 2s, 2p orbitals 99 xi Abstract Quantum materials have emerged as a material of choice for new computing paradigm and information storage. Emergent properties of quantum materials are generally studied under non- equilibrium conditions such as photoexcitation and high electric field. Modeling far from equilibrium condition for quantum materials is a new frontier in molecular modeling. Quantum material’s exotic properties arise from intrinsically quantum phenomena, and thus requiring a quantum treatment of non-equilibrium processes during modeling. This thesis explores the application of quantum molecular dynamics (QMD) under three different non-equilibrium conditions for various materials. First, we perform nonadiabatic quantum molecular dynamics (NAQMD) simulations to study photoexcitation-driven non-thermal amorphization of phase-change materials, GeTe and Sb2Te3. Phase change materials are of high interest for low-power, high-throughput storage devices for next generation neuromorphic computing technologies. Their usage is based on the contrasting properties of their amorphous and crystalline phases, which can be switched in the nanosecond timescale. Here, we show that the energy required for, and the time associated with, the amorphization process can be further reduced by using a non-thermal photoexcitation path. In our NAQMD simulations for GeTe and Sb2Te3, non-thermal pathway was characterized by an instantaneous charge transfer from Te-p orbitals to Ge-p/Sb-p orbitals upon excitation. Time evolution of excited state shows an increasing bonding characteristic for M-M (M = Ge, Sb) and decreasing M-Te interaction, leading to the formation of wrong bonds. These results highlight an energy-efficient and ultrafast amorphization pathway that could be used to enhance the performance of phase change materials based opto-electronic devices. xii Second, we model shock in aramid fibers with quantum and classical molecular dynamics (MD) simulations using multiscale shock theory (MSST). Aramid fibers are widely used as lightweight, shock-resistance, reinforcing materials. However, their intrinsic shock behaviors are not well understood. We perform ab-initio QMD simulations of shock on poly(p-phenylene terephthalamide) (PPTA), which reveals stress release mechanisms based on structural phase transformation (SPT) and generation of planar amorphous region. SPT is triggered by [100] shock- induced coplanarity of phenylene groups and rearrangement of sheet stacking, leading to a previously unknown monoclinic phase. Planar amorphous region is generated by [010] shock- induced scission of hydrogen bonds leading to disruption of polymer sheets, and trans-to-cis conformational change of polymer chains. In contrast to the latter, the former mechanism preserves hydrogen bonding and cohesiveness of polymer chains in the identified novel crystalline phase preserving the strength of PPTA. The interplay between hydrogen bond preserving (SPT) and non- preserving (planar amorphous) shock release mechanisms are critical to understanding the shock performance of aramid fibers. To understand the size effect, we also perform million-atom MD simulation of Kevlar using reactive forcefield. In [100] direction, we observed the structural phase transformation similar to that in QMD. On the contrary, in [010] direction, we found the formation of paracrystalline Kevlar for large systems, which was impossible to study by QMD. Third, we model non-equilibrium shock simulation to understand the detonation process in 2,4,6-triamino-1,3,5-trinitrobenzene (TATB). Detonation processes probed with atomistic details have remained elusive due to highly complex reactions in heterogeneous shock structures. Here, we provide atomistic details of the initial reaction pathways during shock-induced decomposition of TATB crystal using large reactive molecular dynamics simulations based on reactive force fields. Simulation results reveal the existence of three competing intermolecular pathways for the xiii formation of N2. We also observe the formation of large nitrogen- and oxygen-rich carbon aggregates, which delays the release of final reaction products. Finally, we use adiabatic QMD to perform computational synthesis of Nitride ferroelectrics, Al1-xScxN (x = 0.27 – 0.43). We use a Berry-phase approach in conjunction with QMD to study electric field-induced switching of polarization in Nitride ferroelectrics. The computed electric field value for polarization switching in Al0.75Sc0.25N (600 MV/m) is in close agreement with the experimental value (400 MV/m). Furthermore, our simulations correctly predict the lowering of switching field with increasing Sc concentration. QMD also determines the switching pathway and effect of Sc doping in Al1-xScxN. Apart from all the traditional MD techniques, machine learning (ML) has recently grown in popularity. In particular, ML has been used to study the ground state of materials. Here, we use ML to extract both ground and excited states of the hydrogen atom from the Schrödinger equation. The ML technique used to learn Schrödinger equation provides energies and wavefunctions of ground and excited states, which are in excellent agreement with the analytical results. Remarkably, the combined ML and self-consistent methods provide exact wave functions without Gram-Schmidt orthogonalization. With ever-increasing knowledge of materials, we need better tools to perform high- throughput screening of material and model fundamental and exotic properties occurring in the materials with high fidelity. In the present era of quantum materials where non-equilibrium processes are prevalent, this thesis demonstrates a roadmap to study non-equilibrium process using quantum molecular dynamics and machine learning. 1 1. Introduction The ability to compute and store information has been an important aspect of human development. From cave paintings to the discovery of abacus and archival scrolls show the human desire for computation and information storage. Introduction of semiconductors revolutionized the process of computation and information storage. In the recent years, quantum materials promise to be a fundamentally new approach to computation and storages such as quantum and neuromorphic computing. Advances in quantum computing could enable us to solve extremely hard problems such as factoring integers efficiently in minutes.1 Google’s recent paper demonstrated 1.59×109 time speedup by a 53-qubit quantum computer over current state-of-the- art supercomputer.2 Similarly, advancement of neuromorphic computer allows us to mimic learning and adaptation as human brain.3 Quantum materials are defined by solids with exotic properties arising from quantum mechanical nature of their constituent electrons which do not have a classical analogue.4,5 Non- equilibrium conditions such as the application of electromagnetic field and photoexcitation have been used to drive quantum materials far from equilibrium and study their emergent properties.6 Photoexcitation has been used to control magnetic and ferroelectric properties of materials and induce phase change.7–10 One such example is laser-induced heat pulse to write patterns on Skyrmions.11 Electric-field control can also be tailored for properties in demand for quantum materials such as magnetic chiralities in ferroaxial multiferroic RbFe(MoO4)2.12,13 These complex quantum phenomena occurring under non-equilibrium conditions at the sub-picosecond timescale are mainly investigated experimentally. However, molecular dynamics (MD) simulation can be tailored for such unique class of problems. The inherent electronic nature of quantum materials 2 restricts the application of MD to quantum molecular dynamics (QMD). Further, the prevalent occurrence of non-equilibrium processes must be included in the QMD simulation. In this thesis, we explore the modeling of non-equilibrium processes, specifically, electric field, photoexcitation and shock, using QMD simulation. To showcase the applicability of these state of the art non-equilibrium methods, we have performed QMD simulation on various quantum and non-quantum materials: 1) nonadiabatic quantum molecular dynamics (NAQMD) simulation of non-thermal amorphization of phase change materials such as GeTe and Sb2Te3; 2) electric field switching of ferroelectric material, AlxSc(1-x)N; 3) shock response of Kevlar; 4) non-equilibrium dynamics of 2,4,6-Triamino-1,3,5-Trinitrobenzene crystal. Further, we use a new and emerging technique in machine learning (ML) called equation learning to learn the Schrödinger equation of hydrogen atom, where we identify ground and excited states without involving Gram-Schmidt orthogonalization process. This method could potentially speed up costly non-equilibrium quantum-mechanical calculations in QMD on future AI accelerators. In this chapter, we will discuss the brief introduction on these techniques and their application. 1.1. Light driven non-thermal amorphization process in GeTe and Sb2Te3 Photoexcitation is one of the fundamental non-equilibrium processes known to human. In recent year, photoexcitation has been used to control chemical and physical properties of the materials.7,8,10,14,15 Specifically in quantum materials, photoexcitation has shown to control and induce the phase change of materials10 which can be utilized to store information. One such material, Ge2Sb2Te5, transformed the information storage technologies in 20th centuries.16 These materials can be switched from crystalline-amorphous-crystalline within approximately 200 nanosecond using laser induced heating (figure 1a).17,18 Extreme differences between 3 optoelectronic properties of their crystalline and amorphous phase made them ideal for devices. Recently, Ge-Sb-Te alloys gained attention for their application in neuromorphic computing and low power high throughput storage devices.16,19–22 The SET process in Ge2Sb2Te5 alloys is transition from crystalline to amorphous transition which occurs within approximately 10 ns. we propose a faster crystalline to amorphous transition pathway by employing photoexcitation with stronger and short laser pulse. Here, the amorphization of crystal occurs via non-thermal pathway within picosecond compared to previous nanosecond process, as shown in figure 1b. The ultrafast timescale and low energy requirement for SET process via non-thermal amorphization are highly sought after for high transfer rate devices. Figure 1: Scheme Set and Reset process in phase change materials a) via melt-quench process using laser b) via non-thermal photoexcitation. 4 In chapter 3, we study the non-adiabatic quantum molecular dynamics (NAQMD) simulation of non-thermal picosecond amorphization of GeTe and Sb2Te3, main constituent of Ge2Sb2Te5, induced by photoexcitation. Our study explicitly models the instantaneous electronic excitation process, characterizes the evolution of electronic states, their effect on the GeTe and Sb2Te3 bonding, the electron phonon coupling, the diffusion processes, and the resulting loss of long- range order. 1.2. Hydrogen bond is key to shock resilience of Aramid fibers Aramid fibers, such as Kevlar and Twaron, are of great interest due to their outstanding ratio of strength-per-weight and thermal resistance.23,24 There are numerous applications of these fibers that include bullet proof vests and helmets, turbine engine fragment containment structures, cut resistant gloves, etc.25,26 Composition and structure studies of aramid fibers suggest that these materials derive their exceptional properties from crystalline domains of poly (p-phenylene terephthalamide) (PPTA) polymers, highly oriented along the fiber axis.27–32 While aramid fibers display exceptional performance, their development and improvement have been historically empirical and based on trial and error.33,34 Recent investigations on effect of defects on the performance of fibers have been instrumental in the understanding intrinsic deformation and failure mechanism.35–37 However, considering the fact that the main application of these fibers is on shock and impact and protection, there is a startling lack of shock investigations on PPTA and aramid fibers.38–40 Such studies are required for the understanding of the intrinsic shock damage mechanisms of PPTA. Atomistic simulations are key to understanding shock damage microscopically and to enable atomistically-informed continuum simulations.41,42 Recent simulations have highlighted the influence of defects, and the arrangement of chain structures, on the mechanical properties.35– 5 37,43Nevertheless, in order to understand the intrinsic ability of aramid fabrics to withstand shock loading, the response of its most basic constituent, PPTA crystals, needs to be investigated explicitly first. Ideally, the shock studies should describe accurately subtle quantum effects in PPTA and be modeled by first principles methods.44,45 Such first-principles approach to shock based on the Multiscale Shock Technique (MSST) was recently proposed.46 In chapter 4, we will combine MSST shock theory with quantum and classical molecular dynamics to study the shock response of Kevlar. Further we use active learning to identify a new phase in Kevlar. 1.3. Multiple reaction pathways in shocked 2,4,6-triamino- 1,3,5 trinitrobenzene crystal Many organic high explosive (HE) materials are metastable molecular solids composed of C, H, N and O atoms. During detonation, HE materials produce stable small molecules such as N2, H2O and CO2 via complex chemical reactions and release excess energy.47 The time scale of these detonation processes varies widely from nanoseconds to microseconds, depending on the molecular composition and physical properties (e.g., density and crystal structure) of HEs as well as chemical reaction pathways involved in the decomposition processes.48,49 However, these processes are very difficult to understand with atomistic details.50–52 2,4,6-triamino-1,3,5-trinitrobenzene is an extremely insensitive yet moderately powerful polyamino-polynitro arene HE material.53 Its unusual thermal and shock resistances, as well as its chemical stability, make TATB highly important, technologically.54 Due to the interplay of covalent bonds with hydrogen bonds and van der Waals interactions, reactivity of TATB under shock is also a fundamental scientific problem.49,55,56 Understanding of its chemical reaction during 6 detonation process with atomistic details may provide an insight about the chemical stability of TATB and eventually leads to the discovery of safer and more powerful HE materials. In chapter 5, we performed non-equillibribum molecular dynamics simulation to understand the reaction mechanism for TATB decomposition under shock. 1.4. Electric field induced polarization switching in AlxSc(1-x)N Ferroelectricity in materials originates from the spontaneous polarization which is quantum mechanical in nature itself.57 The spontaneous electronic polarization can be switched with application of electric field. This phenomenon has been utilized into RF devices, non-volatile memory, integrated circuit technology, sensors, actuators and neuromorphic computing.58–60 Currently employed traditional ferroelectric, such as HfO2/ZrO2 and perovskites, suffer from very diffusive transition in hysteresis curve.59,61 This leads to a long switching time between two polarization states. Recently discovered Nitride ferroelectrics, Al1-xScxN (x = 0.27 – 0.43), counter such issue and show a sharp square hysteresis curves even in polycrystalline films shown in figure 2 (Taken from Ref 28).62 Further Nitride ferroelectrics, Al1-xScxN (x = 0.27 – 0.43), are highly attractive for device purposes as their spontaneous polarization large (60—100 µC/cm²) compared to tradition ferroelectric such as PZT (20-40) and Fe-HfO2(1-40).62,63 Even though experimentally available, local structure of Al1-xScxN and role of Sc in ferroelectricity is still unknown. Further, the process leading to extremely sharp transition between two state in states are also unknown. In Chapter 6, we will perform quantum molecular dynamics to create Al1- xScxN alloys with different concentration and study the local structure of these alloys. Further, using a Berry-phase approach64–66 with QMD we will model the electric field induced polarization 7 switching in Nitride ferroelectrics. This chapter provide a general framework for computational synthesis and investigation of switching pathways in ferroelectric nitride. Figure 2: hysteresis curves of Al1-xScxN (x = 0.27 – 0.43)62 1.5. Machine learning solution of Schrödinger equation via principle of least action Neural networks have been known for their capabilities as universal function approximators.67 In recent years, significant progress has been made in development of neural networks (NN) which has establish machine learning at the forefront of data driven models. However, in science community, machine learning (ML) remains naïve as many data driven approaches completely ignore prior knowledge within NN models. There are several attempts made to utilize such knowledge to create a better machine learning models.68–70 Raissi et al. 8 proposed a neural network approach to solve ordinary differential equations (ODEs) and partial differential equations (PDEs) with nonlinear form 𝑢 𝑡 +𝒩 [𝑢 ;𝜆 ]=0, 𝑥 ∈Ω,𝑡 ∈[0.𝑇 ] with autodiff71 and physics informed neural networks. Here 𝒩 (𝑢 ;𝜆 ) is a nonlinear operator parameterized by 𝜆 , is a subset of ℝ 𝐷 , and 𝑢 (𝑡 ,𝑥 ) is the latent deep learning solution for the PDE. Raissi et al. have used both continuous and discrete time models to solve various non-linear problems such as the time-dependent Schrödinger equation, Burger’s equation and Navier-Stokes equation. In Chapter 7, We have also used autodiff in conjunction with a self-consistent scheme to solve the Schrodinger equation for the hydrogen atom. This approach provides energies and wavefunctions of ground and excited states that are in excellent agreement with the analytical results. Remarkably, the combined ML and self-consistent methods provide exact wave functions without the use of explicit Gram-Schmidt orthogonalization. 9 2. Method 2.1. Overview of molecular dynamics Molecular dynamics (MD) is an initial value problem.72–74 It follows the physical movements of particles at any given time from their initial positions using Newton’s equation of motion under a given interatomic force field. A N-body molecular dynamics simulation consist of N particles with 3N-element positions, 𝑟⃗ 𝑁 =(𝑥 1 ,𝑦 1 ,𝑧 1 ,𝑥 2 ,𝑦 2 ,𝑧 2 ,……𝑥 𝑁 ,𝑦 𝑁 ,𝑧 𝑁 ) and the velocities, 𝑣⃗ 𝑁 =(𝑣 1𝑥 ,𝑣 1𝑦 ,𝑣 1𝑧 ,𝑣 2𝑥 ,𝑣 2𝑦 ,𝑣 2𝑧 ,……𝑣 𝑁𝑥 ,𝑣 𝑁𝑦 ,𝑣 𝑁𝑧 ) . The trajectories of each atom are predicted by numerically solving the Newton’s equations of motion, where interatomic forces between particles and the associated potential energy are defined by a parameterized force field. 𝑎⃗ 𝑖 (𝑡 )= 𝐹⃗ 𝑖 (𝑡 ) 𝑚 𝑖 (2.1) where 𝑚 𝑖 is the mass of 𝑖 𝑡 ℎ particle. If t is a small increment in time, the velocity of ith particle can be derived from the change of position with respect to time, 𝑣⃗ 𝑖 (𝑡 )= 𝑑 𝑟⃗ 𝑖 𝑑𝑡 = lim ∆𝑡 →∞ 𝑟⃗ 𝑖 (𝑡 +∆𝑡 )−𝑟⃗ 𝑖 (𝑡 ) ∆𝑡 (2.2) The acceleration of the 𝑖 𝑡 ℎ particle can be inferred from the change in velocity with respect to time, 𝑎⃗ 𝑖 (𝑡 )= 𝑑 2 𝑟⃗ 𝑖 𝑑 𝑡 2 = 𝑑 𝑣⃗ 𝑖 𝑑𝑡 = lim ∆𝑡 →∞ 𝑣⃗ 𝑖 (𝑡 +∆𝑡 )−𝑣⃗ 𝑖 (𝑡 ) ∆𝑡 (2.3) 10 We use velocity-Verlet numerical method to integrate the Newton’s equation of motion for each particle. In the velocity-Verlet algorithm, expression for the velocity of each particle in the acceleration form is given by half the time step, 𝑎⃗ 𝑖 (𝑡 )= lim ∆𝑡 →∞ 𝑣⃗ 𝑖 (𝑡 + ∆𝑡 2 )−𝑣⃗ 𝑖 (𝑡 − ∆𝑡 2 ) ∆𝑡 = lim ∆𝑡 →∞ 𝑟⃗ 𝑖 (𝑡 +∆𝑡 )−2𝑟⃗ 𝑖 (𝑡 )+𝑟⃗ 𝑖 (𝑡 −∆𝑡 ) ∆𝑡 2 (2.4) Using equation (2.1) to Eq. (2.4), we can relate the positions and velocities of a system of N particles with given interatomic force fields. Thus, as seen from the spatio-temporal evolution of a physical system represented by that system could be predicted with given initial values. 2.2. Interatomic forcefield To describe the state of each particle in the system, it is essential to compute the force applied on each particle. In molecular dynamics, we consider atomic forces that are derived from the total potential energy, which govern the motion of the atoms. The potential energy can be written as a function of a 3N-dimensional position vector, where N is the total number of atoms present in the system. 𝑉 =𝑉 (𝑟 1 ,…..𝑟 𝑁 )=𝑉 (𝑥 1 ,𝑦 1 ,𝑧 1 ,…..𝑥 𝑖,𝑦 𝑖,𝑧 𝑖,…..,𝑥 𝑁 ,𝑦 𝑁 ,𝑧 𝑁 ) (2.5) Force on the ith atom is represented as 𝐹 𝑖 ⃗⃗⃗ =− 𝜕𝑉 (𝑟 1 ,…..𝑟 𝑁 ) 𝜕 𝑟⃗ 𝑖 = −( 𝜕𝑉 𝜕 𝑥⃗ 𝑖 , 𝜕𝑉 𝜕 𝑦⃗ 𝑖 , 𝜕𝑉 𝜕 𝑧⃗ 𝑖 ) (2.6) In general, the total interatomic potential energy could contain multiple terms corresponding to a specific physical effect. A very commonly used interatomic potential used for 11 inert gases is Lennard-Jones potential. Figure 3 shows a graph between the potential energy and distance in dimensionless units. 𝑉 (𝑟 )=∑𝑉 (𝑟 𝑖𝑗 )= 𝑖 <𝑗 ∑4𝜖 [( 𝜎 𝑟 𝑖𝑗 ) 12 −( 𝜎 𝑟 𝑖𝑗 ) 6 ] 𝑖 <𝑗 (2.7) Figure 3: Dimensionless Lennard-Jones potential for inert gases. 2.3. Reactive forcefield For detonation and shock simulation, we used a variant of reactive force field (ReaxFF)75 called ReaxFF-lg76, which can describe chemical reactions and dynamic properties for large-scale reactive system, with both bonded and non-bonded interactions.76,77 For bonded terms, bond order is applied to take account of the formation and dissociation of the chemical bond. For non-bonded 12 terms, ReaxFF considers shielded van der Waals (Morse Potential) and Coulomb interactions with a finite range as shown in Figure 4.75 Figure 4: Energy terms included in ReaxFF where Ebond is the bond energy, Elp is the lone pair energy, Eover and Eunder are the over- and under- coordination energies, Eval is the valence angle energy, Epen is the penalty energy, Ecoa is the three-body conjugation energy, Etors is the torsion angle energy, Econj is the four-body conjugation energy, Ehbond is the hydrogen bond interaction energy, EvdWaals, ECoulomb and Elg are the van der Waals, Coulomb energies and dispersion energy, respectively. We use ReaxFF-lg, ReaxFF with dispersion correction as shown in equation 2.8. ReaxFF- lg has been tested for several HEM such as RDX, PETN and TATB. ReaxFF-lg provides a relatively better equation of state of TATB than ReaxFF. ReaxFF-lg also provides a better ground state description (such as heat of sublimation, density and cell parameters). 𝐸 𝑅𝑒𝑎𝑥𝐹𝐹 −𝑙𝑔 =𝐸 𝑅𝑒𝑎𝑥𝐹𝐹 + ∑ 𝐶 𝑙𝑔 ,𝑖𝑗 𝑟 𝑖𝑗 6 +𝑅 𝑒 𝑖𝑗 6 𝑁 𝑖<𝑗 (2.8) 13 Where, Reij is equilibrium vdW distance between atom i and j. 𝐶 𝑙𝑔 ,𝑖𝑗 is the dispersion energy correction parameter. 𝑟 𝑖𝑗 6 is distance between two atom i and j. 2.4. Quantum molecular dynamics In Quantum mechanics, Hamiltonian of the N-electrons system can be defined based on the positions of atomic nuclei and electrons. The Hamiltonian of the system can be defined by equation 2.9 𝐻̂ = − ℏ 2 2𝑚 𝑒 ∑∇ 𝑖 2 𝑖 +∑ 𝑍 𝐽 𝑒 2 |𝑟 𝑖 −𝑅 𝐼 | 𝑖,𝐼 + 1 2 ∑ 𝑒 2 |𝑟 𝑖 −𝑟 𝑗 | 𝑖≠𝑗 −∑ ℏ 2 2𝑀 𝐼 ∇ 𝐼 2 𝑖 + 1 2 ∑ 𝑍 𝐼 𝑍 𝐽 𝑒 2 |𝑅 𝐼 −𝑅 𝐽 | 𝐼 ≠𝐽 (2.9) where ℏ is Planck’s constant, MI and ZI are the mass and charge of the Ith nucleus; me and e are the electron mass and charge, respectively. Lower-case indices span the electrons and upper- case indices span the nuclei. Current popularity of quantum molecular dynamics as a predictive tool can be attributed to the development of two methods. 1) development of Pseudopotential 2) development of Exchange correlations. Pseudopotentials (PP) methods are used to describe the electronic wave function, Φ(𝑟 ) . In presence of a periodic potential, motion of electrons should satisfy the Born-von Karman periodic boundary condition as shown in equation 2.10 Φ(r+N j a j )= Φ(𝑟 ) (2.10) Where 𝑗 =1,2,3 and 𝑁 =𝑁 1 𝑁 2 𝑁 3 is the number of primitive unit cell in the crystal; 𝑁 𝑗 is the number of unit cells in the jth direction. Applying the Bloch theorem on periodic potential leads to wavefunctions formed by Bloch wavefunction.78 However, the cost of an electron solution for heavy atoms with Bloch wavefunction would increase exponential. This leads to the pseudopotential method, where atomic nucleus and core electrons are chemically active valence 14 electrons are separated.79 Most pseudopotential currently used i.e. NCPP, USPP satisfies the four general conditions of pseudopotentials. 1) Pseudopotential should contain no node; 2) Normalized atomic-radial pseudo-wave-function with angular momentum is equal to the normalized radial all- electron wavefunction (AE) beyond chosen cutoff; 3) Charges enclosed within chosen cutoff radius must be same and 4) valence all-electron and pseudopotential eigen values must be same.80– 82 In 1994, Projector Augmented Wave (PAW) method was introduced which combined the versatility of linear projector augmented method (LPAW) and PP method. PAW has become state- of-the-art method since its introduction due to its all electron description of the atom at the cost of PP methods.83 Development of exchange-correlations have an extremely rich history. In this chapter, I would like to focus on Density functional theory (DFT), which has become a workhorse for current QMD simulation. Hohenberg-Kohn theorem states that the ground state energy 𝐸 can be expressed as a functional of the ground state electronic density 𝜌 as 𝐸 =𝐸 [𝜌 ]. Thus, the total energy of the system can be approximated in the pseudopotential formalism by equation 2.11. 𝐸 [𝜌 ]=𝑇 [𝜌 ]+ 𝐸 𝑒𝑙 [𝜌 ]+𝐸 𝐻 [𝜌 ]+𝐸 𝑥𝑐 [𝜌 ]+𝐸 𝑙𝑙 [𝜌 ] (2.11) Here, 𝑇 [𝜌 ] is the kinetic energy, 𝐸 𝑒𝑙 [𝜌 ] contains local and non-local term arising from PP formalism developed by Kleinman-Bylander.84 Core and electron interaction can be represented by coulomb interaction. 𝐸 𝑥𝑐 [𝜌 ] is the exchange correlation energy which is defined by equation 2.12 𝐸 𝑥𝑐 [𝜌 ]= 1 2𝜋 ∫ 𝑑 3 𝑟 𝜌 (𝑟 )𝜖 𝑥𝑐 (𝜌 (𝑟 )) 2𝜋 0 (2.12) 15 Figure 5: Workflow schematic for general QMD simulation to obtain Kohn-Sham eigenstates. After input of initial atomic positions {Ri} and atomic velocities {vi}, steps include making an initial guess for the charge density (r); computation of the system Hamiltonian, including kinetic energy T, Hartree potential VH, exchange-correlation potential VXC, local potential VL, and non-local potential VNL; solving the Kohn-Sham equations; computing new charge density ’(r); checking convergence criterion; outputting new atomic positions if criterion is met or mixing ’(r) and (r) to begin to next iteration if not. 𝜖 𝑥𝑐 (𝜌 (𝑟 )) in equation can be approximated as function of electronic density, called local density approximation (LDA).85 Further sophistication has been imposed on approximation of 𝜖 𝑥𝑐 (𝜌 (𝑟 )) , where it becomes function of local density and gradient of density, known as general gradient approximation (GGA). One such function, PBE (Perdew, Burke, Ernzerhof), is commonly 16 used in density functional theory.86 𝐸 𝑙𝑙 [𝜌 ] is the coulomb energy of the all the nuclei as shown in equation 2.13 𝐸 𝑙𝑙 [𝜌 ]=∑ 𝑍 𝑖 𝑍 𝑗 |𝑅 𝑖 −𝑅 𝑗 | 𝑖 <𝑗 (2.13) In the Kohn-Sham (KS) density functional theory (DFT), ground state energy of a many electron interacting system is found by solving an effective N 1 electron noninteracting problem self-consistently.87 Car and Parrinello combined the density functional theory with molecular-dynamics techniques where electronic and atomic part are treated simultaneously while integrating the Newton’s equation of motion.88 Figure 5 shows the workflow for QMD simulation.89 2.5. Non-adiabatic quantum molecular dynamics Quantum molecular dynamics with electronic transition are an important class of problems, which describe the photoexcitation dynamics of electrons and ions. In non-adiabatic quantum molecular dynamics (NAQMD) simulation90–94 based on Tully’s fewest switches surface hopping method,95,96 electrons at each time step occupy the excited eigenstates corresponding to the atomic configurations. The electronic transition between eigenstates is stochastically determined based on the transition probability by a non-adiabatic coupling (NAC). The NAC between excited states is described by a density matrix and its time evolution is calculated using time-dependent density functional theory (TD-DFT).97–101 The ionic motion follows the classical Newtonian mechanics and interatomic forces are computed using Hellman-Feynman theorem. The many-body effect in NAQMD is introduced by describing electronic excited state as a linear combination of electron- hole pairs within the linear-response TDDFT.102–105 In Kohn Sham TDDFT, electronic excited 17 states are expanded in term of ground state Kohn-sham (KS) orbitals as a basic set. However, application of local exchange correlation (xc) functional in KS-DFT leads to several known problems of such as underestimation of band gap, failure to describe charge-transfer excited state. The former problem is corrected by self-interaction correction or GW approximation while the later problem is corrected by range-separated hybrid xc functionals that incorporate long-range exact exchange correlation or by incorporating many-body approach such as Bethe-Salpeter equation (BSE).85,106–110 To reduce the cost of exact exchange functional, NAQMD simulations employ non-self-consistent (NSC) approximation to hybrid functional. 2.5.1. Long-range corrected ground state dynamics Ground state wavefunctions are first obtained self-consistently using GGA functional. Long range exchange correlation through a range-separated hybrid exact exchange (xc) functional is added on self-consistent GGA KS orbitals. The xc energy functional is given by equation 2.14 𝐸 𝑥𝑐 =𝐸 𝑐 ,𝐺𝐺𝐴 +𝐸 𝑥 ,𝐺𝐺𝐴 𝑆𝑅 +𝐸 𝑥 ,𝐺𝐺𝐴 𝐿𝑅 (2.14) where 𝐸 𝑐 ,𝐺𝐺𝐴 is the GGA correlation energy functional, 𝐸 𝑥 ,𝐺𝐺𝐴 𝑆𝑅 is short-range part of GGA exchange energy functional and 𝐸 𝑥 ,𝐺𝐺𝐴 𝐿𝑅 is long-range part HF exchange integral. We obtain the long-range corrected KS energies and orbitals within NSC approximation.111 Long-range interaction is computed using the reciprocal-space formalism by Martyna et al. to avoid interaction with periodically-repeated image.112 2.5.2. Electronic excitations We describe electronic excited states as a linear combination of electron-hole pairs within Casida’s LR-TDDFT using ground state long range corrected KS orbitals as a basis set. In LR- 18 TDDFT, electronic excitation energies are calculated from the pole of an electron-hole pair response function by solving a non-Hermitian eigenvalue problem105 as shown in equation 2.15 ( 𝐴 𝐵 𝐵 ∗ 𝐴 ∗ )( 𝑋 𝐼 𝑌 𝐼 )=𝜔 𝐼 ( 1 0 0 −1 )( 𝑋 𝐼 𝑌 𝐼 ) (2.15) Where 𝜔 𝐼 is the 𝐼 -th excitation energy, with the corresponding eigenvectors 𝑋 𝐼 and 𝑌 𝐼 in equation 2.18. 𝐴 and 𝐵 in equation 2.18 is given by 𝐴 𝑎𝑖𝜎 ,𝑏𝑗𝜏 =𝛿 𝑎 ,𝑏 𝛿 𝑖,𝑗 𝛿 𝜎 ,𝜏 (𝜖 𝑎𝜎 ′ −𝜖 𝑖𝜎 ′ )+𝐾 𝑖𝑎𝑖𝜎 ,𝑏𝑗𝜏 𝐵 𝑎𝑖𝜎 ,𝑏𝑗𝜏 =𝐾 𝑖𝑎𝑖𝜎 ,𝑏𝑗𝜏 Where the indices i, j and a, b are used for occupied and virtual orbitals, respectively. 𝜎 ,𝜏 are spin variables and 𝜖 𝑖𝜎 ′ is the 𝑖 𝑡 ℎ LC-KS orbitals energy with spin 𝜎 .107,108 The coupling matrix 𝐾 𝑖𝑎𝑖𝜎 ,𝑏𝑗𝜏 is given by K ias,jbt = ¢ y is * ¢ y as 1 r é ë ê ¢ y jt ¢ y bt * ù û -d s,t ¢ y is * ¢ y js é ë erf(mr) r ¢ y at ¢ y bt * ù û + dr ò d ¢ r ò ¢ y is * (r) ¢ y as (r) d 2 (E xc,GGA -E x,GGA LR ) dr s (r)dr t ( ¢ r) ¢ y jt ( ¢ r) ¢ y bt * ( ¢ r) (2.16) where E xc,GGA is the xc functional within GGA. The size of the A, B and K matrices in Equations (2.18)-(2.19) is NoNuNoNu, where No and Nu respectively are the numbers of occupied and unoccupied LC-KS orbitals used to represent excited states. According to the assignment ansatz of Casida, the many-body wave function of the I-th excited state is given by 19 F I = iÎ{occupied} å aÎ{unoccupied} å X I,ias +Y I,ias w I ˆ c as + ˆ c is F 0 s å (2.17) where F 0 is the Slater determinant of the occupied LC-KS orbitals, and ˆ c ks + and ˆ c ks are the creation and annihilation operators acting on the i-th LC-KS orbital of spin . 2.5.3. Force in molecular dynamics For an excited electronic state, we use NSC method similar to the Harris-Foulkes approach to evaluate accurate excited-state forces at a moderate computational cost. Here, the ground state NSC Harris-Foulkes energy ENSC is defined as E NSC = f n y n ˆ H KS y n y n y n n å -E dc (2.18) where Edc represents double-counting terms related to the Hartree and xc energies: E dc =- 1 2 dr ò r(r)V H (r)+E xc r [ ] - dr ò r(r)V xc (r) (2.19) ENSC coincides with the KS energy EKS for self-consistent charge density, which is generally determined by the error bounding the difference between in(r) and out(r). The force Fk acting on the k-th atom is derived from ENSC in Equation (2.18) as F k =- ¶ ¶R k E NSC + z l z m R l -R m l<m å æ è ç ç ö ø ÷ ÷ =F k Hellmann-Feynman +F k NSC (2.20) where Rk and zk are the position and the valence, respectively, of the k-th atom. In Equation (2.20), FkNSC is the NSC force, which appears when the self-consistency is not satisfied, and is given by F k NSC = drd ¢ r dr( ¢ r) r- ¢ r - ¶ ¶R k r in (r) ì í î ü ý þ òò + drdr(r) - ¶ ¶R k V xc (r) ì í î ü ý þ ò (2.21) 20 where δ(r) = out(r) − in(r). For self-consistent charge density, FkNSC =0 as δ (r) = 0. We can rewrite I-th excited state ΦI, equation (2.17), in the LR-TDDFT method as F I = C I,ias ias å ˆ c as + ˆ c is F 0 º C I,ias ias å F ias (2.22) where Φ i aσ is the Slater determinant obtained by replacing the i-th orbital by the a-th orbital in Φ0. From this simplification, excitation energy ωI is formally written as w I = F I ˆ H N F I (2.23) where ˆ H N is the Hamiltonian of the N-electron system. Atomic forces in the excited state are usually calculated by finite difference method. To reduce the computational cost of calculating excited-state forces, we first employ a diagonal approximation defined as F k =- ¶ ¶R k w I + z l z m R l -R m l<m å æ è ç ç ö ø ÷ ÷ ~- ¶ ¶R k C I,ias ias å 2 F ias ˆ H N F ias + z l z m R l -R m l<m å æ è ç ç ö ø ÷ ÷ (2.24) In the spirit of the NSC force, we evaluate the first term of Equation (2.24) using the LC-KS orbitals ¢ y n and their occupation numbers f n ias in Φiaσ as − 𝜕 𝜕 𝑅 𝑘 ⟨Φ 𝑖𝑎𝜎 |𝐻̂ 𝑁 |Φ 𝑖𝑎𝜎 ⟩= 𝜕 𝜕 𝑅 𝑘 {∑𝑓 𝑛 𝑖𝑎𝜎 𝑛 ⟨𝜓 𝑛 ′ |𝐻̂ 𝑘𝑠 |𝜓 𝑛 ′ ⟩ ⟨𝜓 𝑛 ′ |𝜓 𝑛 ′ ⟩ −𝐸 𝑑𝑐 } (2.25) Here, we have used 𝜌 𝑖𝑎𝜎 ′ (𝑟 )= ∑𝑓 𝑛 𝑖𝑎𝜎 𝑛 |𝜓 𝑛 ′ | 2 (2.26) 21 as out(r), and the electron density of the ground state as in(r) to evaluate the value of F k NSC,ias .113 2.2.6 Nonadiabatic electron-ion dynamics In NAQMD simulations with electronic transitions using Tully’s FSSH method95,96 along with the LC-KS representation of TDDFT, we calculate the time evolution of the density matrix with a fixed atomic configuration between MD steps. The elements of the density matrix determine the switching probability between the adiabatic states that are expressed by a single Slater determinant. The equations for the density matrix are solved using the eigen energies of the single- electron LC-KS Hamiltonian and the NAC terms, which are obtained from the current (excited) adiabatic state. These equations are derived by expanding the electronic state Y(t) at time t in terms of the electronic excited states F J (R(t)) in LR-TDDFT corresponding to the atomic configuration R(t) at time t: Y(t) = C J (I) (t) F J (R(t)) J å C J (I) (0)=d I,J (2.27) The time evolution of the expansion coefficients C J (I) (t) is governed by d dt C J (I) (t)=- C k (I) (t) iw K d JK +D JK ( ) k å (2.28) where the NAC elements are defined as D JK = F J ¶ ¶t F K (2.29) 22 From a pair of consecutive time steps in an adiabatic QMD simulation, the NAC elements provide phonon-assisted electronic transitions. 2.6. Multiscale shock Theory (MSST) The MSST technique allows simulation of shock waves with significantly lower cost. The shock wave is simulated by uniform deformation of crystalline unit cells subject to macroscopic conservation laws of mass, momentum and energy. In multiscale shock theory, propagation of shock wave is modeled by 1D Euler’s equation for compressible gas flow. Euler equation for 1D compressible gas law ignores the thermal transport. 1D compressible gas equation corresponding to mass momentum and energy conservation is given in equation 2.30. 114 𝑑𝜌 𝑑𝑡 +𝜌 ( 𝜕𝑢 𝜕𝑥 )=0 2.30(a) 𝑑𝑢 𝑑𝑡 +𝜈 ̅( 𝜕𝑃 𝜕𝑥 )=0 2.30(b) 𝑑 𝑒 ̃ 𝑑𝑡 + 𝑃 ( 𝜕 𝜈 ̅ 𝜕𝑡 )=0 2.30(c) Where 𝜌 is mass density, 𝑢 = 𝑑𝑟 𝑑𝑡 is local materials velocity, 𝑟 is position of material, 𝜈 ̅ is specific volume. 𝑃 is pressure and 𝑒 ̃ is energy per unit mass. As the shock passes through the system, we follow the position of atoms in the frame of reference of shock wave. Thus, position of ith particle (xi) in new frame of reference given by (xi -vst). where 𝑣 𝑠 is shock speed. Using equation 2.30 and position of particle in new frame of reference we can obtain three new equation of state. 23 𝑢 −𝑢 0 =(𝑢 −𝑣 𝑠 )(1− 𝜌 𝜌 0 ) 2.31(a) 𝑃 −𝑃 0 =(𝑢 𝑜 −𝑣 𝑠 ) 2 𝜌 0 (1− 𝜌 𝜌 0 ) 2.32(b) 𝑒 ̃−𝑒 ̃0 =𝑃 0 ( 1 𝜌 0 − 1 𝜌 )+ (𝑢 𝑜 −𝑣 𝑠 ) 2 2 (1− 𝜌 𝜌 0 ) 2 2.32(c) Equation 2.32 are called Rankine-Hugoniot conditions for conservation of mass, momentum and energy. Equation 2.14(b) is Rayleigh line. In the laboratory frame of reference material is at rest, we can assume 𝑢 0 =0. Using equation 2.32 and replacing 𝜌 =1/𝑉 and 𝜌 0 =1/𝑉 0 , lagrangian for MSST equation can be derived. 𝐿 = 1 2 ∑𝑚 𝑖𝑟 ̇ 2 𝑖 −𝜙 ({𝑟 𝑖 })+ 1 2 𝑄 𝑀 𝑉 ̇2 + 1 2 𝑀 𝑣 𝑠 (1− 𝑣 𝑠 𝑣 0 ) 2 +𝑃 0 (𝑉 −𝑉 0 ) (2.33) Where Φ is the potential energy; 𝑀 =∑𝑚 𝑖 𝒊 is the total mass of the system. mi is the mass of the ith atom. Q has units of 𝑚𝑎𝑠 𝑠 4 /𝐿𝑒𝑛𝑔𝑡 ℎ 4 . 𝑉 0 is the initial volume of the system and 𝑉 is the volume of system at time t. In this work, we used the omni-directional multiscale shock technique (OD-MSST), implemented by Shimamura et al., which is a generalization of the MSST technique.46,56 Equation of motion of atoms is modified by dynamically evolving the volume of crystalline unit cell, while constraining the stress to the Rayleigh line and energy to the Hugoniot equation. The extended Lagrangian of the system is given in equation 2.34 𝐿 = 1 2 ∑𝑚 𝑖 (ℎ𝑞 𝑖̇) 𝑡 (ℎ𝑞 𝑖̇)−𝜙 ({ℎ𝑞 𝑖 },ℎ)+ 𝑄 2𝑀 𝑉 2 + 1 2 𝑀 𝑣 𝑠 2 (1− 𝑉 𝑉 0 ) 2 − 𝑃 0 (𝑉 −𝑉 0 ) 𝑖 (2.34) 24 qi is a column vector whose components are the ith atom’s scaled coordinates in the range [0, 1], and vS is the speed of the shock wave. The real coordinate and the velocity of the ith atom are given by hqi and ℎ𝑞 𝑖̇ , respectively, where h = (L1 L2 L3) is a matrix consisting of the computational cell lattice vectors Lk (k = 1, 2, 3). V = det (h) is the volume of the computational cell. P0 and V0 = det(h0) are the pressure and volume of the unshocked state, respectively, where h0 corresponds to h in the unshocked state. In Eq. (2.34), a dot denotes time derivative. MSST can be implemented with both quantum molecular dynamics and classical molecular dynamics. 2.7. Finite electric field in quantum molecular dynamics In presence of electric field with periodic boundary condition, the total energy of the system can be defined by variational energy functional as shown in equation 2.35.64 𝐸 𝐸 [{𝜓 𝑖 }]= 𝐸 0 [{𝜓 𝑖 }]−Ω𝐸 𝑃 [{𝜓 𝑖 }] (2.35) where, 𝐸 0 [{𝜓 𝑖 }] is the ground state energy of the system, Ω is volume of system and 𝐸 is applied electric field. 𝑃⃗⃗ [{𝜓 𝑖 }] is polarization calculated using the Berry phase method as shown in equation 2.36.115,116 𝑃 [{𝜓 𝑖 }]= − L 𝜋 Im(ln (𝑑𝑒𝑡 [𝑀 ])) (2.36a) ℳ=⟨𝜓 𝑖|𝑒 −2𝑖𝜋𝑥 /𝐿 |𝜓 𝑗 ⟩ (2.36b) The Hamiltonian is then constructed as follows: 𝐻 |𝜓 𝑖 ⟩=𝐻 0 |𝜓 𝑖 ⟩−𝑉 𝑐𝑒𝑙𝑙 ℇ ⃗⃗⃗⃗ ∙ 𝛿 𝑃⃗⃗ [{𝜓 𝑖 } 𝛿 ⟨𝜓 𝑖 | where: (2.37) 25 𝛿 𝑃⃗⃗ [{𝜓 𝑖 }] 𝛿 ⟨𝜓 𝑖 | =−𝑖 𝑒 2𝐿 2 𝜋 (∑|𝜓 𝑗 ⟩𝑀 𝑗𝑖 −1 −∑|𝜓 𝑗 ⟩(𝑀 𝑗𝑖 ∗ ) −1 ) 𝑗 𝑗 where M is defined in equation 2.36. To calculate this gradient, we have used 𝑑 𝑑𝑥 Log Det 𝐴 = Tr 𝑑𝐴 𝑑𝑥 𝐴 𝑗𝑖 −1 and Im(𝑧 )= −𝑖 2 (𝑧 −𝑧 ∗ ) . In the crystalline system approach this gradient becomes: 𝛿 𝑃⃗⃗ [{𝜓 𝑛 𝒌 𝒋 }] 𝛿 ⟨𝜓 𝑛 𝒌 𝒋 | =− 𝑖𝑒 2Δ𝑘 (∑|𝑢 𝑚 𝑘 𝑗 ⟩𝑆 𝑚𝑛 −1 (𝒌 𝒋,𝒌 𝒋 +𝟏 )−∑|𝑢 𝑚 𝑘 𝑗 −1 ⟩𝑆 𝑚𝑛 −1 (𝒌 𝒋,𝒌 𝒋−𝟏 ) 𝑚 𝑚 𝑆 𝑛𝑚 (𝒌 𝒋,𝒌 𝒋 +𝟏 )=⟨𝑢 𝒌 𝒋 𝑛 |𝑢 𝒌 𝒋 +𝟏 𝑚 ⟩ Here, 𝑢 𝒌 𝒋 𝑛 is cell periodic part of 𝜓 𝑛𝑘 . With these gradients, a general optimization technique can be used to find the ground state in the presence of a finite electric field. (2.38) (2.39) 26 3. Light Driven Non-thermal Amorphization in Sb2Te3 and GeTe 3.1. Introduction With ever-increasing demand for faster electronics and physical and manufacturing constraints of current technologies, there is a growing interest in new materials that enable hardware acceleration and are suitable for new computing paradigms.16,22,117 Phase change materials (PCMs) are one of the most promising materials to fulfill the role of both high throughput non-volatile memory devices16,19,20 and brain inspired neuromorphic computing.21,22 Devices based on PCMs leverage the contrasting optoelectronic properties between their amorphous and crystalline phases.118–120 Ge-Sb-Te alloys, an archetypical phase change material, have been widely investigated both experimentally and theoretically due to their rapid crystalline-amorphous transition.16,119 Specifically, Ge2Sb2Te5 can be switched from crystalline to amorphous in short time span of 10 ns by applying short laser or electric pulses. On the other hand, the Ge2Sb2Te5 crystallization process requires an order of magnitude longer time.16,22 Other Ge-Sb-Te alloys has been explored such as GeSb2Te4, GeSb3Te7 and Ge3Sb2Te6 along the quasibinary line of GeTe and Sb2Te3.121 Speed of amorphous to crystalline transition (SET) and crystalline to amorphous transition (RESET) is controlled by changing the GeTe/ Sb2Te3 ratio. With decreasing Sb2Te3, SET temperature increases and with increasing Sb2Te3 SET temperature decreases. GeTe and Sb2Te3 themselves are also phase change materials. Sb2Te3 has shown a much faster phase change cycle and lower power consumption compared to Ge2Sb2Te5, due to its growth dominated crystallization and lower melting point.17,122 On the contrary, GeTe crystallization and amorphization occurs both in hundreds of nanoseconds. However, timespan of RESET (crystalline to amorphous transition) can 27 be further reduced by picosecond timescale using photoexcitation process by stronger and short laser pulses, i.e. femtosecond laser. While being technologically important process, the underlying mechanism and change in chemical bonding for non-thermal amorphization of all these materials are still not understood. Photoexcitation amorphization is a common phenomenon where materials undergo amorphization as a result of electronic excitation followed by electron-phonon coupling, gradual heat up of the lattice and melting followed by quenching.16,118 This method is employed in rewritable CD, DVD and Blu-ray technologies. By employing strong femtosecond laser pulses, high electronic excitations can be achieved leading to the instability of materials crystalline structure.18,123 That opens a path for direct non-thermal amorphization, i.e. loss of long range order without involving a melt-quench process.124–126 Time-resolved optical measurements and theoretical studies report ultrafast sub-picosecond amorphization induced by femtosecond lasers in different materials.119,120,123,127 The ultrafast time-scale and low energy involved in the non- thermal amorphization process are highly sought after for energy efficient and high transfer rate devices.[citation] Previous theoretical studies estimate that about 11% valence electron excitation is required to induce lattice instability in covalently bonded silicon and gallium arsenide.18,128 The crystalline structure of GeTe, Sb2Te3, as well as other Ge-Sb-Te alloys, is stabilized by a network of resonant bonds. It is still to be determined, the exact electronic and ionic processes leading to the destabilization of GeTe and Sb2Te3 crystalline structure and resulting non-thermal amorphization by femtosecond laser photoexcitation. Here, we present a non-adiabatic quantum molecular dynamics (NAQMD) simulation of non-thermal picosecond amorphization of GeTe and Sb2Te3 induced by photoexcitation. Our study explicitly models the instantaneous electronic excitation process, characterize the evolution of electronic states, their effect on the GeTe and 28 Sb2Te3 bonding, the electron phonon coupling, the diffusion processes, and resulting loss of long- range order. 3.2. Non-thermal amorphization of Sb2Te3 3.2.1. Simulation setup To systematically study the effects of photoexcitation in Sb2Te3 and the possibility of non- thermal amorphization, we perform NAQMD simulations on a supercell consisting of 180 atoms with n = 2.6, 5.2, 10.3, and 15% valence electrons in excited states, corresponding to a carrier density of 0.45, 0.91, 1.82, and 2.44 (1022 cm-3), respectively. Figure 6a illustrates the initial Sb2Te3 crystalline supercell. The crystalline structure is initially thermalized at 300 K using Nose- Hoover thermostat for 600 fs shown in Figure 6b. After initial thermalization, we excite the valence electrons and simulate the evolution of excited states at constant temperature of 500 K within canonical ensemble (NVT). The rationale for this simulation setup is to highlight the effect of photoexcitation on the structure stability of Sb2Te3, assuming that the slow buildup of heat from local photoexcitation is readily dissipated throughout a large experimental sample. Figures 6c-6f show the structure of Sb2Te3 after photoexciting n = 2.6, 5.2, 10.3 and 15.0 % valence electrons 2.6% (Figure 6c) do not affect the crystalline order of the structure within the time evolution considered. In contrast, for n = 10.3% excitation, Figure 6e indicates loss of long-range order, which also applies for n = 15% (Figure 6f). Considering that the melting temperature of Sb2Te3 is 890 K129, the result suggests an amorphization of the structure via a path not involving the usual melting-quenching, i.e. a non-thermal amorphization. For n = 5.2%, Figure 6d indicates a local rearrangement of the structure while most of the structure remains crystalline. 3.2.2. Analysis of crystalline to amorphous phase change To obtain further insights into the excited dynamics leading to the contrasting results illustrated in Figure 6c-f, we calculate the mean square displacement (MSD) of all aforementioned 29 cases. The MSD results for Sb atoms are shown in Figure 7. As a reference point, we performed a simulation with n = 0% as shown by the black line in Figure 7, which indicates no atomic diffusion. Results for n = 2.6%, shown by the cyan line in Figure 7, closely follows the reference curve for n = 0%, indicating no change in structure. In contrast, the green curve in Figure 7, corresponding to the simulation with n = 5.2% excitation, shows a residual MSD of 0.6 Å2, which is clearly higher than that of n = 0% and n = 2.6%. A visual analysis of Figure 6d suggests that the higher residual MSD is a result of picosecond induced localized disorder in the structure. Results from simulations with higher excitation levels, n = 10.3% and 15%, shown by the orange and red curves, indicate a distinct behavior characterized by a sharp increase in the MSD within 2 ps timescale followed by steady increase of the MSD up to 4Å2 within 4 ps. That indicates that electronic excitation at these levels induces large atomic diffusion. A visual analysis of Figure 6e, for n = 10.3% excitation, shows that the induced atomic diffusion leads to widespread bond breaking and loss of the crystalline arrangement and long-range order. Figure 6: Sb2Te3 simulation supercell structures (a) initial crystalline structure of Sb2Te3, (b) structure at 300K (c-f) structure of Sb2Te3 after 4.22 ps excited state dynamics with n = 2.6, 5.2, 10.3 and 15.0 % excited valence electrons, respectively. Yellow and purple spheres correspond to antimony and tellurium atoms, respectively 30 Figure 7: Mean square displacement of Sb as a function of time. Cyan, green, orange and red curve corresponds to excitation of n=2.6%, n=5.2%, n= 10.3% and n=15% respectively at temperature 500 K. Black curve corresponds to the adiabatic MD at 500 K. n is percentage of excited electron out of the total electrons. To further identify the change in the structure after excited state dynamics, we perform radial distribution function analysis. Figure 8a shows the radial distribution function of systems for after 4.22 ps of time evolution. The radial distribution function’s first peaks are located at 3.02Å and 4.32 Å. The former indicates the Sb-Te bonds. The latter indicates the Sb-Sb and Te- Te nearest distance. While these two peaks are well defined for n = 0%, 2.6% and 5.2%, a merged single broad peak is shown for the n = 10.3% and 15% suggesting a change in the short-range order. To understand the change in short range we evaluate the fraction of wrong bonds formed during the evolution of excited state as shown in Figure 8b. For, n = 0% we can observe that the fraction of wrong bonds remains 0%. Beside small fluctuations the same applies to n = 2.6%. As we increase the excitation, the number of wrong bond formed increases. For n = 5.2%, the fraction of wrong bonds remains below 10% suggesting a limited change in short range order. This is reflected in the increase in value of g(r) between first and second peak. For higher excitation (n = 31 10.3% and 15%), the number of wrongs bonds increases dramatically up to 40% and the above suggesting a loss of short-range order. This observation is directly correlated with the merging of the first and second peaks in g(r). Figure 8: Analysis of the structural evolution (a) total radial distribution function as a function of number of excitons (b) evolution of the number of wrong bonds for different excitons (c) evolution of simulated diffraction pattern for n= 10.3% excitation. A loss of long-range order is suggested in Figure 8a. The absence of distinct 3rd, 4th and 5th peak of the crystalline structure shown by the g(r) curves for n = 10.3% and 15% is an indication of loss of long-range order. Calculated electron diffraction pattern, Figure 8c, for the system as a function of time further supports the loss of long-range order and formation of an amorphous system. Figure 8c shows the calculated electron diffraction pattern as a function of time for n = 10.3%. We can observe a strong peak for the [015] plane at 28.3˚ angle as well as lower intensity 32 peaks for [1010], [110] and [205] planes for the initial crystalline system (at time t = 0 ps). At t = 0.59 ps after photoexcitation, the [015] peak is still well defined, while peak intensity for the remaining peaks of the crystalline structure declined significantly. Finally, at t = 4.22 ps, we observe the complete disappearance of the [015] peak suggesting a complete loss of long-range order and formation of an amorphous state. 3.2.3. Evolution of excited state In order to understand the loss of crystalline order and in particular the loss of short-range order, we turn to the changes and evolution of electronic states induced by the photoexcitation process. We evaluate the charge transfer between atomic species by calculating the average Mulliken charge per atomic species during the excited state dynamics. The results for n = 2.6%, 5.2%, and 10.3% are shown in Figure 9. In all cases, the plots show a sharp charge transfer from Te to Sb at time t = 0 ps, the instantaneous photoexcitation time. Negative/positive time in the plots indicate dynamics prior/post photoexcitation. At n = 2.6% excitation, the results in Figure 9a show a nearly identical Mulliken charge for both species after excitation. The evolution of the excited states indicate that the ground state is mostly recovered in a short period of 4 ps. Given that the structural result from Figure 7 and Figure 8 suggest a crystalline state, we expect a full recovery of the ground state for both Sb and Te species in sufficiently longer relaxation. The results for n = 5.2% excitation, shown in Figure 9b, indicates an almost reversed average charge after photoexcitation compared to their ground state. Similar to the n = 2.6% case, we observe a swift partial recovery of the Sb and Te ground state charge within 4 ps of evolution. In contrast, Figure 9c shows a different behavior for n = 10.3%. At such high excitation level, the average charge reverse sign and increase in magnitude indicating a strong charge transfer from Te to Sb. In a short timespan of 4 ps, a strong reversal in charge transfer occurs, even though it indicates the final 33 steady state is distinct from the ground state. Similar effect is observed for n = 15% as shown in Figure 4d. Figure 9: Average Mulliken charge per component as a function of time. (a) n=2.6%, (b) n = 5.2%, (c) n = 10.3% excitation (c) n = 15% excitation. To understand the effect of charge transfer from Te to Sb on the crystalline order, we perform Mulliken bond overlap analysis as a function of time as shown in Figure 10. Figure 10 a- c show the histogram for Mulliken bond overlap for n = 2.6% photoexcitation for Sb-Sb, Sb-Te, and Te-Te bond pairs, respectively. Before photoexcitation, at t = -0.5ps, the reference bond overlap histograms are very well defined, corresponding to the different bonds present in the Sb2Te3 crystal structure. These are represented by the black curves in Figure 10 a-c. Figure 10a indicates a peak at -0.15 for Sb-Sb bond pair, corresponding to antibonding interaction between Sb-Sb atom types. For Sb-Te bond type (Figure 10b), we observe a histogram peak at 0.4, corresponding to a large Sb-Te bonding interaction. For Te-Te bond type, we observe a large peak at -0.07 and a small peak at 0.12. The antibonding peaks at -0.07 corresponds to the intralayer 34 repulsion between Te atoms while the bonding peak at 0.12 corresponds to the interlayer interaction between Te atoms. Such interlayer bonding interaction between Te atoms is associated with the layered structure of Sb2Te3. The same black reference curves are shown in Figure 10 for all other excitation. After excitation at time t = 0 ps, we do not observe any significant changes in Sb-Te bond overlapping histogram for all excitation levels, shown by red curves. In contrast, Sb-Sb histogram peaks for n = 2.6% and n = 10.3% are shifted towards increased antibonding interaction, i.e., negative shift by 0.03 to 0.05. The negative shift of the peak location for Sb-Sb bond types increases with increasing carrier density. In contrast to changes in Sb-Sb and Sb-Te bond types, the Te-Te peak intensity in the negative region decreases and increases in positive region, indicating an enlarged bonding interaction between Te atoms. For n = 2.6%, as the system evolve, we observe that all the aforementioned changes in excited states revert back to the ground state as shown by the green, cyan, and blue curves in Figure 10 a-c. That is consistent with the preservation of the Sb2Te3 crystalline structure. However, for n = 10.3%, as the system evolves, the Sb-Te bond overlap curve becomes broad and expands into the antibonding region shown by the green, cyan, and blue curves in Figure 10h. That corresponds to the scission of Sb-Te bonds. Concurrently, the two peaks displayed by Te-Te at time t = 0 ps, merged into a broad distribution indicating no clear bonding and antibonding interaction (Figure 5i). The results for Sb-Sb, displayed in Figure 10g, indicates that the antibonding interaction is still prevalent although residual bonding interaction is present. Similar process is observed for n = 5.2 and 15% (Figure 5d-f and5j-l). The combination of Sb-Te bond scission and Sb-Sb/Te-Te bond formation is consistent with the amorphization of the structure as illustrated in Figure 6d. 35 Figure 10: Evolution of Mulliken bond overlap analysis as a function of time for a) Sb-Sb, b) Sb-Te, and c) Te-Te bond pairs for n = 2.6%. Similarly, d-f) correspond to Sb-Sb, Sb-Te and Te-Te bond pairs for n = 5.2%, g-i) correspond to Sb-Sb, Sb-Te and Te-Te bond pairs for n = 10.3%, and j-l) correspond to Sb-Sb, Sb-Te and Te-Te bond pairs for n = 15%, respectively. Black 36 color curve corresponds to the t = -0.5 ps (before excitation). Red, green, cyan, and blue colors correspond to times t = 0 ps , 0.5 ps, 1.0 ps, and 2.5 ps after excitation. 3.3. Non-thermal amrophization of GeTe 3.3.1. Simulation setup To systematically study the effects of photoexcitation in GeTe crystal, we perform NAQMD simulations on a supercell consisting of 144 atoms with n = 3.3, 5, 7.5, and 10% excited electrons, corresponding to a carrier density of 6.08, 9.07, 13.6, and 18.1 (1021 cm-3), respectively. Figure 11a illustrates the initial GeTe crystalline supercell. The crystalline structure is initially thermalized at 300 K using Nose-Hoover thermostat for 600 fs shown in Figure 11b.130,131 After initial thermalization, we excite the valence electrons and simulate the evolution of excited states at constant temperature of 800 K within canonical ensemble (NVT).130,131 The rationale for this simulation setup is to highlight the effect of photoexcitation on the structure stability of GeTe, assuming that the slow buildup of heat from local photoexcitation is readily dissipated throughout a large experimental sample. Figure 11c-e show the structure of GeTe after photoexciting n = 3.3, 5, 7.5 and 10.0 % valence electrons. For n = 3.3% (Figure 11c), we observe a local rearrangement of the structure while most of the structure remains crystalline. In contrast, for n = 5% excitation, Figure 11d, indicates loss of long-range order, which also applies for n = 7.5% (Figure 11e). Considering that the melting temperature of GeTe is 1000 K,129 the result suggests an amorphization of the structure via a path not involving the usual melting-quenching, i.e. a non- thermal amorphization. 37 Figure 11: GeTe simulation supercell structures (a) initial crystalline structure of Sb2Te3, (b) structure at 300K (c-f) structure of GeTe after 4.22 ps excited state dynamics with n = 2.6, 5.2, 10.3 and 15.0 % excited valence electrons, respectively. Orange and blue sphere represent germanium (Ge) and tellurium (Te) atoms, respectively 3.3.2. Amorphous to crystalline transformation To obtain further insights into the excited dynamics leading to the contrasting results illustrated in Figure 11c-e, we calculate the mean square displacement (MSD) of all aforementioned cases. The MSD results for Ge atoms are shown in Figure 12. As a reference point, we perform a simulation with n = 0% as shown by the black line in Figure 12 at 800K temperature, which indicates no atomic diffusion. Results for n = 3.3%, shown by the cyan line in Figure 12, closely follows the reference curve (n = 0%) with small residual MSD, indicating a small rearrangement. In contrast, the green curve in Figure 12, corresponding to n = 5% excitation, shows a residual MSD of 1.2 Å2, which is clearly higher than that of n = 0% and n = 3.3%. A visual analysis of Figure 11d suggests that the higher residual MSD is a result of disorder in the structure. 38 However, it is important to note that almost flat line MSD after 1ps suggest an amorphous phase formation. Results from simulations with higher excitation levels, n = 7.5% and 10%, shown by the orange and red curves, indicate a distinct behavior characterized by a sharp increase in the MSD within 0.6 ps timescale followed by steady increase of the MSD up to 4Å2 within 2.5 ps. That indicates that electronic excitation at these levels induces large atomic diffusion. A visual analysis of Figure 11e, for n = 7.5% excitation, shows that the induced atomic diffusion leads to widespread bond breaking and loss of the crystalline arrangement and long-range order. Figure 12: Mean square displacement of Ge atom type as a function of time. Cyan, green, orange and red curve corresponds to excitation of n=3.3%, n=5%, n= 7.5% and n=10% respectively at temperature 800 K. Black curve corresponds to the adiabatic MD at 800 K. n is percentage of excited electron out of the total electrons. 39 Figure 13: Analysis of the structural evolution (a) total radial distribution function as a function of number of excitons (b) evolution of the number of wrong bonds for different excitons (c) evolution of simulated diffraction pattern for n= 10.3% excitation. To further identify the change in the structure after excited state dynamics, we perform radial distribution function analysis. Figure 13a shows the radial distribution function of systems for after 2.5 ps of time evolution. The radial distribution function’s first peaks are located at 2.9 Å and 4.1 Å. The former indicates the Ge-Te bonds and the latter indicates the Ge-Ge and Te-Te nearest distance. While these two peaks are well defined for n = 0%, the separation between Ge- Ge and Te-Te peaks closes creating a shoulder for the n = 5% and 7.5% and 10%. This suggests a change in the short-range order. To understand the change in short range we evaluate the fraction of wrong bonds formed during the evolution of excited state as shown in Figure 13b. For, n = 0% we can observe that the fraction of wrong bonds remains 0%. For n = 3.3%. the fraction of wrong 40 bonds remains below 10% suggesting a limited change in short range order. This is reflected in the increase in value of g(r) where we can observe a separate first and second peak. For higher excitation (n = 7.5, and 10%), the number of wrongs bonds increase dramatically up to 40% and above, suggesting a loss of short-range order. This observation is directly correlated with the merging of the first and second peaks in g(r). For n = 5%. the fraction of wrong bonds rapidly increases to 12% and remains almost constant suggesting an amorphization in first 0.5 ps after which small local rearrangement occurs due to thermal fluctuation. A loss of long-range order is suggested in Figure 13a. Calculated electron diffraction pattern, Figure 13c, for the system as a function of time further supports the loss of long-range order and formation of an amorphous system. Figure 13c shows the calculated electron diffraction pattern as a function of time for n = 5%. We can observe a strong peak for the [200] plane at as well as lower intensity peaks for [111] and [220] planes for the initial crystalline system (at time t = 0 ps). At t = 1.6 ps after photoexcitation, the [200] peak is still well defined, while peak intensity for the remaining peaks of the crystalline structure declined significantly. Finally, at t = 2.5 ps, we observe the significant decrease in intensity of the [200] peak suggesting a complete loss of long- range order and formation of an amorphous state. 3.3.3. Electronic and ionic transformation during excited state dynamics In order to understand the loss of crystalline order and in particular the loss of short-range order, we again turn to the changes and evolution of electronic states induced by the photoexcitation process. We evaluate the charge transfer between atomic species by calculating the average Mulliken charge per atomic species during the excited state dynamics. The results for n = 3.3%, and 10% are shown in Figure 14. In all cases, the plots show a sharp charge transfer from Te to Ge at time t = 0 ps, the instantaneous photoexcitation time. Negative/positive times in 41 the plots indicate dynamics prior/post photoexcitation. At n = 3.3% excitation, the results in Figure 14a show a identical and opposite Mulliken charge for both species after excitation. The evolution of the excited states indicate that the ground state is mostly recovered in a short period of 4 ps. Given that the structural result from Figure 12 and Figure 13 suggest almost crystalline state, we expect a recovery of the crystalline state with small local deformation in sufficiently longer relaxation. In contrast, Figure 14b shows a different behavior for n = 10%. At such high excitation level, the average charge reverse sign. In a short time, span of 4 ps, a strong reversal in charge transfer occurs again, even though it indicates the final steady state is distinct from the ground state. Figure 14: Average mulliken charge per component as a function of time. (a) n=3.3%, (b) n = 10% valence band excitation excitation. To understand the effect of charge transfer from Te to Ge on the crystalline order, we perform Mulliken bond overlap analysis as a function of time as shown in Figure 15. Figure 15 a- c shows the histogram for Mulliken bond overlap for n = 5% photoexcitation for Ge-Ge, Ge-Te, and Te-Te bond pairs, respectively. Before photoexcitation, at t = -0.5ps, the reference bond overlap histograms are very well defined, corresponding to the different bonds present in the GeTe 42 crystal structure. These are represented by the cyan curves in Figure 15 a-c. Figure 15a indicates a peak at -0.21 for Ge-Ge bond pair, corresponding to antibonding interaction between Ge-Ge atom types. For Ge-Te bond type (Figure 15b), we observe a histogram peak at 0.4 and a small shoulder at 0.05, corresponding to a large Ge-Te bonding interaction. For Te-Te bond type, we observe a large peak at -0.1 corresponding to antibonding interaction. Figure 15: Evolution of Mulliken bond overlap analysis as a function of time for a) Ge- Ge, b) Ge-Te, and c) Te-Te bond pairs for n = 5%. Cyan color curve corresponds to the t = -0.1 ps (before excitation). orange, green curves correspond to times t = 0.1 ps , 2.0 ps after excitation. 5d) average mulliken bond overlap for Ge-Ge (red curve), Te-Te(blue curve) and Ge-Te(black curve) as a function of time. After excitation at time t = 0.1 ps, we observe significant changes in Ge-Ge bond overlapping histogram shown by orange curve, which shows a broad peak in bonding region at 0.4 shown by Figure 15a. In contrast, Ge-Te bond overlap curve becomes broad and expands into the 43 antibonding region shown by the orange in Figure 15b. This indicates to a scission of Ge-Te bonds. We observe a decay in antibonding peak of Te-Te in Figure 15c. However, there is no shift in position of peak suggesting antibonding interaction shielded by Ge. As the excited state evolves, Ge-Ge bonding interaction remains prevalent. However, small Ge-Ge antibonding interaction start to appear suggesting formation of Ge-Ge bond. For t= 2.0 ps, we observe further decay in Te-Te antibonding peak. Overall, we can observe that as the time progress Ge-Te bonding character declines while Ge-Ge bonding character increases (Figure 15d) suggesting Ge-Te bond scission and Ge-Ge bond formation. This process is consistent with the amorphization of Ge-Te as illustrated in Figure 11c. 3.4. Conclusion In this chapter, we have found a non-thermal amorphization pathway for GeTe and Sb2Te3 via photoexcitation. The ground state partial density of state analysis of both GeTe and Sb2Te3 seems similar. For both materials, valence band maximum is composed of Te p orbitals while excited states are made of M p orbitals (M = Ge, Sb). After excitation, we observe a charge transfer from Te to M. However, the excited state dynamics differ significantly. We can observe from Sb- Te bonding peak intensity and position does not change upon excitation (Figure 10) and Ge-Te peak decline significantly upon excitation (Figure 15). Further, we observe a significant shift in Ge-Ge antibonding peak to bonding, no such change occurs in Sb-Sb. On the contrary, we observe change in Te-Te peaks for Sb2Te3 in Mulliken bond overlap analysis which is not present in GeTe. These changes suggest completely two different pathways for the amorphization of GeTe and Sb2Te3. 44 4. Shock Simulation of Aramid Fibers 4.1. Introduction Aramid fibers, such as Kevlar and Twaron, display outstanding ratio of strength-per-weight and thermal resistance.23,24 These properties have been highly utilized towards application in bullet-proof vests and helmets, turbine engine fragment containment structures, cut resistant gloves, etc.25,26Composition and structure studies of aramid fibers suggest that these materials derive their exceptional properties from crystalline domains of poly (p-phenylene terephthalamide) (PPTA) polymers, highly oriented along the fiber axis.27–29,31,32 Figure 16a illustrates the aramid fiber structure. While aramid fibers display exceptional performance, their development and improvement have been historically empirical and based on trial and error.33,34 Recent investigations on effect of defects on the performance of fibers have been instrumental in the understanding intrinsic deformation and failure mechanism.35–37 However, considering the fact that the main application of these fibers is on shock and impact and protection, there is a startling lack of shock investigations on PPTA and aramid fibers.38,132 Such studies are required for the understanding of the intrinsic shock damage mechanisms of PPTA. Here we use QMD to understand the intrinsic chemical changes in PPTA under shock loading. Further, we use reactive molecular dynamics (RMD) to perform large-scale long-time simulation. RMD simulation help us understand the effect of simulation box size on result, long time rearrangement of the polymer chains under shock loading, formation in polycrystalline Kevlar. Chapter 4.2 will discuss the quantum molecular dynamics simulation. Chapter 4.3 will discuss the reactive forcefield fitting for Kevlar and Chapter 4.4 will discuss the application of newly parameterized forcefield on RMD shock simulation of Kevlar. 45 Figure 16: Structure of aramid fibers and PPTA (p-phenylene terephthalamide) crystal. (a) Aramid fibers microstructure is composed of elongated, highly crystalline fibrils where PPTA polymer chains are aligned along the fibril length forming stacked PPTA sheets. (b) PPTA polymer chains are held together by in-plane hydrogen bonding (orange colored dotted lines) to create a PPTA sheet. Axes a, b, and c correspond to [100], [010], and [001] crystallographic directions. (c) PPTA crystal oriented along b, highlighting the PPTA sheet stacking structure. (d) 46 PPTA crystal oriented along the polymer chain (fibril length). Grey, white, red, and blue spheres correspond to carbon, hydrogen, oxygen, and nitrogen atoms, respectively. 4.2. MSST simulation with quantum molecular dynamics The models used for shock simulations along the crystallographic directions [100] and [010] are composed of 4×2×1 and 2×4×1 repeating PPTA crystalline units27 in an orthorhombic simulation box with size 31.5 × 10.4 × 12.9 Å3 and 15.7 × 20.7 ×12.9 Å3, respectively. Both simulated systems contain 448 atoms. Initially, the system is relaxed using a Quasi-Newton method (BFGS algorithm) at 0 GPa pressure. After relaxation, the multiscale shock technique (MSST) is used to simulate the system under shock condition. Each simulation is performed for 20000 steps. Time-step for all the simulation is 0.2 a.u. 4.2.1. Simulation details In this work, we investigate the shock response of PPTA using the multiscale shock technique (MSST) in the framework of density functional theory.114 MSST simulations are performed perpendicular to PPTA chains along the [100] and [010] crystallographic directions, at several particle velocities in the range 0.24 - 4 km/s. The PPTA crystal structure is highly anisotropic along the low index crystallographic directions [100], [010], and [001]. Covalent bonding between monomers form long polymer chains along the [001] direction (shown in Figure 16b). In-plane hydrogen bonding along the [010] direction forms polymer sheets (hydrogen bonding is indicated by the dotted lines in Figure 16b). In the [100] direction, PPTA sheets are stacked in ABAB arrangement and interact by weaker van der Waals interactions (shown in Figure 16c and Figure 16d). Therefore, distinct types of bonding exist in each direction of the crystal, i.e. covalent along the chain axis, hydrogen bonds within sheets, and van der Waals between different sheets. This bond hierarchy leads to highly anisotropic mechanical properties, i.e. high tensile strength and 47 modulus along chain axis contrasts with the corresponding compressive properties and with the mechanical properties along other directions.24,133,134 This anisotropy of mechanical properties requires a careful consideration in the design of the fibers. The manufacturing process of commercial fibers, such as Kevlar and Twaron, take advantage of the PPTA crystal structure by producing highly aligned polymer chains along the fiber axis. Figure 17: Shock Hugoniot along [100] (a) and [010] directions (b), respectively. PPTA structure in elastic (c), transformation (d), and cross-linking (e) regimes for shocks along [100] direction and figure (f-h) along [010] direction, respectively. Chains undergoing structural or trans-cis transformation are highlighted. Figure (i-k) Analysis of structural changes in different regimes by the evolution of dihedral angle between phenylene and amide groups in the elastic 48 regime, fraction of hydrogen bonds preserved in the plastic/transformation regimes, and fraction of sp2 preserved in the cross-linking regime. 4.2.2. Hugoniot curve The shock Hugoniot of the PPTA crystal along the [100] and [010] directions are presented in Figure 17a and Figure 17b. Corresponding pressure plots are displayed in Figure 18, a and b. Discontinuous changes in shock Hugoniot curves correspond to transitions between different shock response regimes. For shock direction [100] we observe discontinuity at Up = 2.07 km/s and Up = 2.66 km/s. The former corresponds to a PPTA structural phase transformation (SPT) and the latter to cross-linking of the polymer chains. Along [010], a first discontinuity at Up = 0.4 km/s indicates a transition from elastic-to-plastic regime and the second discontinuity at Up = 3.24 km/s is related to polymer chain cross-linking. 49 Figure 18: Shock-induced pressure in PPTA. (a) Induced pressure during shock along the [100] and (b) [010] directions. Red, blue, and green symbols indicate elastic, inelastic (transformation/plastic), and cross-linking regimes. 4.2.3. Structural transformation For weak shocks, the PPTA crystal response shows a highly anisotropic elastic compression regime along different directions. For shock along [100], we observe high volume reduction leading to compaction of the polymer sheets and conformational changes in the PPTA polymer chains. The volume is reduced up to 29% at 2.07 km/s. The distance between PPTA sheets is reduced from 4.6 Å to 3.9 Å. For elastic shock along [010], we observe a relatively low volume reduction up to 6% at 0.4 km/s before plasticity is triggered. In unshocked PPTA, the two phenylene groups are arranged at opposite dihedral angle with the amide group as shown in Figure 16b (label for dihedral is given in supporting information Figure 19). The dihedral angle for unshocked PPTA is 26º. In the elastic regime, the dihedral angle evolution as a function of time is shown in Figure 17i for shock along both [100] (in blue) and [010] (in red) directions. In the case of [100] direction, the dihedral angle reduces from 26º to 16º at 2.07 km/s. In contrast, the dihedral angle along [010] remains unchanged. 50 Figure 19: Dihedral angle between phenylene group and amide. We computed dihedral angle between the phenylene and amide groups as a function of time to observe the conformational changes during elastic shock For strong shocks, we observe cross-linking of the PPTA polymer chains with the formation of new bonds between chains, which drastically changes its chemical structure. Similar phenomenon was reported for ethane clusters under impact resulting in polymerization.135 We have quantified this process by monitoring the integrity of the PPTA backbone, where all carbon atoms have sp2 hybridization. Figure 16k shows the fraction of sp2 carbon present in the system as a function of time in both cases. For [100] direction, we observe a drastic drop in the fraction of sp2 carbon at Up = 3.42 km/s. A similar behavior is observed for [010] direction at Up = 3.7 km/s. structural transformation under mild shock The PPTA displays an intriguing shock response in the intermediate shock regime 1.44 < Up < 3.24 km/s along [010]. We observe a plastic response characterized by: i) a massive scission of hydrogen bonds; ii) rotation of whole polymer chains; iii) independent rotation of phenylene and amide groups; iv) change in the conformation from trans-to-cis; and v) misalignment of polymer chains. The trans-to-cis conformational change, shown in Figure 20, a and b, has a particular deleterious effect on the PPTA crystal, which distorts and stresses the polymer 51 backbone, affecting its ability to regenerate hydrogen bonding with neighboring molecules in the sheets. This shock release mechanism is akin to the shock induced trans-to-gauche isomerization reported for self-assembled hydrocarbon monolayers. All the listed plastic mechanisms act to destabilize the polymer sheets as shown in Figure 16g.136 Figure 16j shows up to an 80% drop in the number of hydrogen bonds. This plastic response effectively destroys the long-range order of the crystalline alignment of the polymer chains leading to a planar amorphous phase, where polymer chains are weakly bonded by van der Waals interactions and contain both trans and cis conformations. The disordered arrangement of the polymer chains is quantified by the radial distribution function, g(r), of their centers of mass position shown in Figure 20c. Sharp peaks in g(r) are lost in the planar amorphous phase indicating lack of crystalline arrangement. Figure 20: trans (a) and cis (b) conformations of a PPTA polymer chain. PPTA crystals have 100% trans conformation. 20% of polymer chains transform to cis during shock along the [010] direction, leading to a planar amorphous phase. (c) Radial distribution function of polymer 52 chains center of mass in the (001) plane for crystalline (C), planar amorphous (A), compressed crystalline (c-C) and transformed (T) systems. The insets show the corresponding center of mass configurations. In sharp contrast to the shock response along [010], the intermediate shock regime along [100], between Up = 2.07 km/s and 2.66 km/s, leads to a pressure induced SPT that preserves the integrity of the PPTA polymer sheets and the trans conformation of the polymer chains. SPT is a common shock release mechanism in strongly bonded materials.137 Results indicate that for some polymer chains the phenylene group attached to the carbonyl end of the amide group reverses its orientation with respect to the amide group making both the phenylene groups coplanar. This conformational change induces displacement of polymer sheets in both [010] and [001] directions and the stacking rearrangement propagates to neighboring polymer sheets. Since these simulations are computationally costly, the shock response is simulated up to 4.88 ps. During this short time period the identified transformation is incomplete for the particle velocities considered here. Please see Figure 16d for a partially transformed system at Up = 2.2 km/s. To uniquely identify the structure of the transformed phase, we perform global energy minimization by sampling the phase space using the Bayesian optimization method. Initial random search showed that the potential energy surface of the system is too flat, making it challenging to find a global minimum. Figure 21A shows the energy landscape as a function polymer chains’ displacement, which shows multiple minima. Furthermore, computation on each structure is extremely costly. Thus, we employ Bayesian optimization technique to efficiently find the optimal solution.138,139 53 4.2.4. Global structure search using active learning Bayesian optimization (BO) optimizes a black box objective function which is typically expensive to compute due to the amount of time required, monetary cost or an opportunity cost. This class of machine learning problem is suitable for solving problems in continuous domain of maximum 20 dimensions. BO builds a surrogate model for the objective function and quantifies the uncertainty in the surrogate using gaussian process regression. We perform Bayesian optimization of smaller system consisting 224 atoms. Further, we only allow the chain as a whole to move in [010] and [001] direction while preserving the Figure 21: (a) Energy landscape of the compressed system as a function of polymer chains’ reduced displacement in [010] and [001] directions. Reduced displacements are calculated using Gödel encoding (see Supporting information). (b) Predicted new phase of PPTA under shock along the [100] direction. All the phenylene groups in polymer chains are coplanar in this PPTA phase. 54 (c) Unshocked PPTA monoclinic primitive unit cell. (d) Predicted new phase PPTA monoclinic primitive unit cell. hydrogen-bond. These constraints allow us to reduce the phase space from 448 (224 atoms in [010] and [001] direction) to 8-dimensional space. The feasible set space in the present problem is define as a set of 8 tuple values {𝑦 𝑖,𝑧 𝑖 } 𝑤 ℎ𝑒𝑟𝑒 𝑖 =1,2,3,4. The polymer chains are allowed to displace ±0.5𝑦 𝑖,±0.5𝑧 𝑖 in the y and z directions. In Bayesian optimization method, acquisition function selects the next point to be search based on trade-off between exploration and exploitation. We use Expected improvement (EI) as acquisition function. After each iteration of BO, a new structure is created and minimized based on the suggestion. Each computed structure is added to the dataset and maximum EI is computed for next structure. Pseudo-code for Bayesian optimization process is given in algorithm 1. ===================================================================== 1 Initial dataset:{(𝑦 𝑖,𝑧 𝑖 ),𝑓 𝑖},𝑖 =1…𝑛 2 Build a gaussian process prior on 𝑓 . 3. Bayesian optimization: 𝑓𝑜𝑟 𝑗 =1 𝑡𝑜 𝑛 𝑚𝑎𝑥 a. Obtain next (𝑦 𝑗 ,𝑧 𝑗 ) by optimizing acquisition function EI over GP as (𝑦 𝑗 ,𝑧 𝑗 )= 𝑎𝑟𝑔𝑚𝑖𝑛 (𝑦 𝑗 ,𝑧 𝑗 ) 𝐸𝐼 ((𝑦 𝑗 ,𝑧 𝑗 )|{(𝑦 𝑖,𝑧 𝑖 ),𝑓 𝑖 } 𝑖=1..𝑛 ) b. Obtain the ground truth 𝑓 𝑗 by running one round of black box DFT package evaluation for obtaining the potential energy value. c. Obtain a new augmented set {(𝑦 𝑖,𝑧 𝑖 ),𝑓 𝑖 } 𝑖=1..𝑛 +1 n = n+1 d. Update the gaussian process prior on 𝑓 . 4. Expected Improvement acquisition function EI a. (𝑦 𝑗 ,𝑧 𝑗 )= 𝑎𝑟𝑔𝑚𝑖𝑛 𝐸𝐼 𝑛 ({(𝑦 𝑖,𝑧 𝑖 ),𝑓 𝑖 } 𝑖=1..𝑛 ) Here 𝑓 𝑖 is a gaussian process prior obtained by fitting the n data points. b. EI = {(𝜇 − 𝑓 𝑚𝑖𝑛 − 𝜉 )Φ(𝑧 )+ 𝜎𝜙 (𝑧 ), 𝑤 ℎ𝑒𝑟𝑒 Φ and 𝜙 𝑎𝑟𝑒 𝑐𝑢𝑚𝑢𝑙𝑎𝑡𝑖𝑣𝑒 𝑎𝑛𝑑 𝑝𝑟𝑜𝑏𝑎𝑏𝑖𝑙𝑖𝑡𝑦 𝑜𝑓 𝑠𝑡𝑎𝑛𝑑𝑎𝑟𝑑 𝑛𝑜𝑟𝑚𝑎𝑙 𝑧 𝑧 = 𝜇 −𝑓 𝑚𝑖𝑛 −𝜉 𝜎 55 =============================================================== Algorithm 1: Bayesian optimization to compute global minimum We performed 40 Bayesian optimization runs to sample the input space and found several possible structures. The Gaussian process regression is fitted initially to ground truth values and augmented each step to include the new data point and newest evaluated values. The model converges to optimum in nearly 40 iterations as shown in Figure 22 and successfully exploit as well as explore the function space. Figure 22: Bayesian optimization to find global minimum. Our input space is an 8-dimensional space which is hard to depict. We employ a modified Gödel encoding to map 8-dimensional space to 2-dimension space. Δ𝑦 = 𝑙𝑜𝑔 2 (3 𝑦 1 +5 𝑦 2 +7 𝑦 3 +9 𝑦 4 ) , here 0≤𝑦 𝑖 ≤1 Δ𝑧 = 𝑙𝑜𝑔 2 (3 𝑧 1 +5 𝑧 2 +7 𝑧 3 +9 𝑧 4 ) , here 0≤𝑧 𝑖 ≤1 Where, displacement in each polymer chain in y and z direction is dented by y={y1, y2, y3, y4,} z={z1, z2, z3, z4}. 56 4.2.5. Analysis of predicted structure We evaluate the energy landscape of PPTA in the compressed state as a function of individual polymer chains displacement in [010] and [001] directions. The calculated energy landscape is shown in Figure 21a. Energy values are plotted as a function of reduced displacements, Δ𝑦 and Δ𝑧 , corresponding to [010] and [001] directions, calculated using Gödel encoding. In order to investigate the effects the exchange-correlation functional on total energies of both the optimized phases, we optimized both phases in PBE, SCAN functional.86,140 To include the effect of long range vdW interaction, we also optimized the geometries with PBE-DFT-D3 and SCAN-rVV10 functional.141–143 All 4 methods predict that new transformed phase has lower energy compared to the compressed PPTA crystal. Figure 23 shows the difference in energy per unit cell between compressed transformed system and compressed crystal PPTA as a function of compression in [100] direction. Both GGA and meta-GGA functional shows that transformed phase have lower energy. 57 Figure 23: Energy difference of transformed state: We computed energy difference between transformed PPTA state and compressed PPTA state at different uniaxial compression in [100] direction. Transformed state shows lower energy for all used functionals. The predicted structure displays the previously described rotation of one phenylene group. The Phenylene group attached to the carbonyl (oxygen) end of amides flip their orientation and become coplanar to the phenylene group attached to the amine (nitrogen) end of amides, (geometry shown in Figure 21b). Symmetry analysis of the predicted structure indicates that it belongs to the same monoclinic space group as the original PPTA crystal P21/c. Figure 21c and Figure 21d show the primitive cells of the original and transformed systems. The primitive cell of the transformed system shows 90˚ rotation of the hydrogen bonded PPTA sheets and coplanarity of the phenylene groups. Given the similarity in space group, the g(r) of the transformed structure, shown in Figure 20c, contrasts with the g(r) of the crystalline structure. We observe a new peak at 2.8Å indicating the highly compact stacking. The predicted anisotropic shock response of PPTA along the independent directions perpendicular to the polymer chain axis offers a framework to understand the outstanding behavior 58 of aramid fibers and its failure mechanism. Our results show that the SPT induced by shock in [100] preserves the in-plane hydrogen bonding which is linked to the observed strength of PPTA. In contrast, the complete breakdown of hydrogen bonding for shock in [010] shifts the failure mechanism from chain rupture to chain pullout, as identified in polyethylene fibers.36 This process is expected to significantly reduce the tensile strength of PPTA. The effect of hydrogen bonding in the PPTA crystal is similar to -sheet crystal domains in silk.144 Considering our results, each randomly oriented fibril may present widely different shock response, i.e. a [100] oriented fibril will outperform [010] oriented fibril under compressive shock. Since the latter arguably will lose its strength in the plastic regime which start at modest particle velocities of Up ~ 0.4 km/s. This loss of strength linked to hydrogen bond disruption may explain the fibrillation failure process identified as a major failure mode in Kevlar fibers.36 Fibrillation of aramid fibers under impact is a predominant failure mechanism that is not well understood. 4.3. Reactive forcefield modeling of Kevlar Ab-initio molecule dynamics situation of Kevlar shows describe the shock response of PPTA crystals. QMD simulation accurately depicts the bond breaking, bond formation and other chemical transformation occurring the system. However, QMD simulation are extremely small in size and do not include the effect of size. Further QMD simulation cannot be run for a long time. Thus, we need to perform classical molecular dynamics simulation to understand the large scale, long time simulation using ReaxFF forcefield.75 ReaxFF is a classical forcefield which is parameterized using first principle. Further, ReaxFF can describe the bond breaking and formation making it ideal for shock simulations. 59 4.3.1. Development of ReaxFF-lg Forcefield We have refined the ReaxFF-lg, a variant of ReaxFF with dispersion correction for non- bonded interaction,76 reactive force field parameters to accurately describe bond ruptures and bulk behaviors of the PPTA crystal by re-optimizing the bond and van der Waals (vdW) terms for C-N and H-N interactions. Our refined ReaxFF-lg forcefield has improved full bond dissociation curves for the C-N and H-N bonds in the amide linkage along the PPTA chain. Furthermore, the refined ReaxFF-lg describes thermally stable PPTA crystal at 10 K and 300 K. ReaxFF force field parameters for C/H/O/N elements were previously developed by Strachan et al. to describe chemical processes of high energy (HE) materials such as RDX.145 Liu et al. extended the ReaxFF formalism by including a long-range-dispersion term (Elg) to improve intermolecular interactions for HE materials (ReaxFF-lg).76 Using the ReaxFF-lg reactive force field parameters, Mercer et al. performed RMD simulations for PPTA crystals, and studied crystallite failure via covalent and hydrogen bond ruptures in the PPTA structures.37 However, the published ReaxFF-lg force field parameters were not optimized against any quantum mechanics (QM) data of PPTA system. For this reason, we investigated two aspects of ReaxFF-lg force field: 1) C-N and H-N bond dissociation energies for a single Kevlar chain, 2) Intermolecular interactions between single Kevlar chains (i.e., bulk characteristics of the Kevlar crystal structure). Unfortunately, we found that C-N and H-N bonds in the PPTA chain were not consistent with our density functional theory (DFT) calculations, which should be calibrated to study bond ruptures of the PPTA for impact and/or shock simulations. Therefore, we decided to refine bond energy (i.e., two-body interactions) and van der Waals (vdW) terms for C-N and H-N interactions 60 in the ReaxFF-lg reactive force field parameters. Those parameters were re-optimized using the one-parameter parabolic search method towards minimizing the total error of the training set, defined as follows: 𝐸𝑟𝑟𝑜𝑟 = ∑ [ (𝑥 𝑖 ,𝑄𝑀 −𝑥 𝑖 ,𝑅𝑒𝑎𝑥𝐹𝐹 ) 𝜎 ] 2 𝑛 𝑖=1 (4.1) where xi,QM and xi,ReaxFF are the QM and ReaxFF values, respectively, and 𝜎 is the assigned weight that is defined in the training set. In the training set, we included QM data of full bond dissociation curves for C-N and C-H bonds in the amide linkage along the PPTA chain using B3LYP/6-31G** using Qchem software.146 Please note that the optimization process above was implemented for the ReaxFF stand-alone code, and thus, we excluded Elg terms in ReaxFF-lg while re-optimizing the bond and vdW terms. After that, we reset Elg terms to the re-optimized parameters. We then tested our refined ReaxFF-lg reactive force field for C-N and C-H bond strengths in the single PPTA chain, and bulk characteristics of the PPTA crystal structure. The results of the force field refinement are discussed in the below. 4.3.2. Bond dissociations energy The central C-N bond in the amide linkage was identified as the weakest bond.37 Thus, it is essential to reasonably describe bond strengths in the amide linkage along the PPTA chains. To 61 Figure 24: Results of the ReaxFF refinement for the Kevlar systems: (a) the central C-N bond in the amide linkage before (right) and after (left) refining; (b) the C-N bond between the amide linkage and aromatic ring in the single Kevlar chain before (right) and after (left) refining; (c) the central H-N bond in the single Kevlar chain before (left) and after (right) refining. Note that energy contributions of the long-range-correction terms (Elg) were not included in the full bond dissociation curves above. 62 correctly capture full bond dissociation curves for C-N bonds as well as H-N bond in the amide linkage, we included the data above in the training set. As shown in Figure 24, the initial ReaxFF- lg force field parameters (left) were not able to reproduce full bond dissociation curves for the C- N and H-N bonds by our DFT calculations, whereas the refined ReaxFF force field parameters (right) quantitatively mimic the bond dissociation energies. We then re-included Elg terms in the refined parameters to reset to the ReaxFF-lg formalism and evaluated bond strengths for the C-N and H-N bonds above as shown in figure 25. We further confirmed that the refined ReaxFF-lg force field parameters still reproduce the bond strengths for the C-N and H-N bonds correctly. Figure 25: Bond dissociation energies described by the refined ReaxFF force field, DFT and ReaxFF-lg: (a) C-N single bond; (b) C-N double bond; (c) H-N single bond. Note that after the refinement, the ReaxFF force field still reproduces the C-N single/double bonds and H-N single bond, which were included in the previous training set. 63 It should be noted that the initial ReaxFF force field parameters for C-N and C-H interactions were primarily trained against DFT calculations using NH2CH3 (a single C-N bond), NHCH2 (a double C-N bond), and NH3 (a single N-H bond) clusters.145 Therefore, we also included full dissociation curves for these clusters in the training set to ensure that the refined force field description still keeps the previous agreement between ReaxFF and DFT for the single/double C- N bonds and single N-H bond unchanged. As a result, we were able to describe these interactions reasonably well, Figure 25. 4.3.3. Energy volume curve We validate the refined ReaxFF-lg force field parameters against DFT for Energy volume curve of Kevlar. For validation, we have taken the geometries starting geometries from QMD MSST simulation (chapter 4.2.1). The initial systems are compressed and expanded in [100] and [010] direction up to 30%. Final energies are computed by optimizing both geometries with DFT and ReaxFF-lg. Energy-volume curve for planar compression of Kevlar in [100] and [010] direction is shown in Figure 26. From Figure 26, we can observe that ReaxFF-lg shows an error of approximately 1 kcal/mol per monomer units for both compressed and expanded system compared to DFT. Overall our refined ReaxFF-lg parameters describe the system with sufficiently good accuracy compared to DFT values. Figure 26: Energy volume curve of Kevlar after planar compression in [100] and [010] 64 Figure 27: RMD simulations of thermal equilibrations for the Kevlar crystal structure: (a) at 10 K for 20 ps; (b) at 300 K for 125 ps. During the thermal process, the Kevlar structure kept the AB arrangements stable in the xy plane. In addition, no significant distortion or bending was observed in the amide linkage (the xz plane) for up to 300 K. 4.3.4. Characterization of PPTA crystal Using the refined ReaxFF-lg force field parameters, we evaluated the equilibrium lattice parameters of the PPTA crystal at 300 K. To obtain a thermally equilibrated PPTA crystal, we 65 used 3×3×3 unit cells (in the [100], [010], and [001] directions, respectively) and performed NPT- RMD simulations. The simulation schedule consists of the following order: In Table 1 below, we show the comparison of the equilibrium lattice parameters between the refined ReaxFF-lg description and values from references. The results suggest that the refined ReaxFF-lg description properly bulk characteristic, consistent with the reference and our DFT calculations. Table 1: Calculated equilibrium lattice constants by DFT, ReaxFF-lg, and refined ReaxFF-lg (this work), compared with experimental values. Equilibrium lattice constants (Å) a b c Expt. (Northolt) 27 7.87 5.18 12.90 Expt. (Rao et al.)32 7.90 5.20 12.92 DFT (this work) 7.49 5.12 13.13 ReaxFF-lg (Liu et al.)76 7.64 4.94 13.33 Refined ReaxFF-lg (this work) 7.61 4.98 13.49 In addition, thermal stability of the PPTA crystal were evaluated for the refined ReaxFF- lg force field. As shown in figure 3a-b, the PPTA crystal structure kept the AB arrangements stable in the xy plane at 10 K and 300 K(right), and no significant distortion or bending effect was observed in the amide linkage at 10 K and 300 K (left). 66 4.4. MSST simulation with Reactive molecule dynamics We performed MSST simulation of Kevlar with refined ReaxFF-lg forcefield to observe the effect of large-scale long-time simulation. We validated our results against QMD simulations. Our MSSST simulation with reactive forcefield shows close match with QMD simulation for small system. Increasing the size introduce the effect of deformation and reorganization of the polymer chain, which is not captured in QMD simulation. 4.4.1. Simulation setup To study the finite size effect of polymer chains, we performed MSST simulations by varying the length of polymer chains in [001] direction from 52.52 Å to 472.68Å while keeping the polymer chains length fixed in [100] and [010] direction. We found elastic to plastic transition point saturated at 149.9 Å. The models used for shock simulations along the crystallographic directions [100] and [010] are composed of 4×2×1 and 2×4×1 repeating PPTA crystalline units To study the effect of large-scale long-time shock simulation, we created models by replicating PPTA crystalline units along the crystallographic directions [100] and [010]. The models are composed of 42×14×28 and 8×77×24 repeating PPTA crystalline units in an orthorhombic simulation box with size 323.89×78.33×378.33 Å3 and 60.79× 394.07×324.05 Å3, respectively. Simulated systems contain 1058400 and 860160 atoms, respectively. Initially, the system is relaxed using a conjugate gradient algorithm. System was further relaxed at 0 GPa pressure using NPT ensemble at 300K. After relaxation, the multiscale shock technique (MSST) is used to simulate the system under shock conditions. Each MSST simulation is performed for 50 ps. 67 4.4.2. Finite size effect of simulation cell Molecular dynamics simulation, performed under a finite simulation cell with periodic boundary condition, generally suffer from size effect. Structural and dynamic properties obtained from MD simulation results depend significantly on size of simulation box. For example, increasing self-diffusivity of liquid atom with increasing system size. Other example is simulation cell size dependent thermodynamics properties.147–151 In MSST shock simulation, effect of system size was observed during transformation of fused silica/quartz to shtishovite.152 Figure 28: shock speed in the Kevlar after immediately Hugoniot shock limit as function of increasing size. To isolate the effect of simulation cell size, we performed MSST simulation of Kevlar in [010] direction with increasing simulation cell length until the Hugoniot elastic limit saturate to a 68 specific particle velocity. We plot shock speed during particle transformation (speed at which system undergoes elastic to plastic transformation after Hugoniot elastic limit) as a function of length of the system in Figure 28. We can observe that DFT (𝑈 𝑠 = 6.7 km/s for elastic-plastic transformation) and ReaxFF (𝑈 𝑠 = 7.0 km/s for elastic-plastic transformation) values are extremely close for small systems. As we increase the simulation cell size to 149 Å, the elastic limit is lowered to 5.8 km/s. Further increase in size generally does not lower the shock speed. Figure 29: Shock Hugoniot along [100] (a) and [010] directions (b), respectively for RMD simulation. 4.4.3. Shock Hugoniot MSST simulation are performed perpendicular to PPTA chains along [100] and [010] crystallographic direction at velocity range of 0.05-3.2 km/s. The shock response of PPTA is shown in Figure 29. We observe at discontinuity at Up = 1.6 km/s and Up = 2.5 km/s for RMD simulation compared to DFT discontinuity at Up = 2.07 km/s and Up = 2.66 km/s in [100] direction. Second discontinuity corresponding to the crosslinking of polymers remains almost unchanged for ReaxFF-lg and DFT within bound or errors. However, we observe a significant changed between first discontinuity corresponding to structural phase transformation. For larger system with 69 ReaxFF-lg, structural transformation occurs at Up = 1.6 km/s compared to DFT (Up = 2.07 km/s). This decrease can be attributed to finite size effect. For [010] direction, finite size effect is much more pronounced where elastic to plastic transformation occurs at 0.05 km/s compared to Up = 0.4 km/s. Second discontinuity appear at Up = 2.9 km/s for ReaxFF-lg compared to Up = 3.04 km/s for DFT. For weak shocks, the PPTA crystal response shows a very similar response to our DFT simulation results. For shock along [100], we observe high volume reduction leading to compaction of the polymer sheets and conformational changes in the PPTA polymer chains. For elastic shock along [010], we observe almost no volume reduction before plasticity is triggered. For strong shocks, we observe cross-linking of the PPTA polymer chains which is very similar to DFT simulation. 4.4.4. Paracrystalline transformation under mild shock The PPTA displays an intriguing shock response in the intermediate shock regime 1.2 < Up < 2.9 km/s along [010]. We observe a short time planar amorphization followed by long time polymer chain rearrangement. Short time plastic response characterized by: i) a massive scission of hydrogen bonds; ii) rotation of whole polymer chains; iii) independent rotation of phenylene and amide groups; iv) change in the conformation from trans-to-cis; and v) misalignment of polymer chains which is very similar to DFT simulation. Long-time simulation shows the complete rotation of polymer chain between 0˚-90˚ and formation of crystalline domain. The boundaries between these domains are not well defined and could span up to 10 rotated polymer chains suggesting the formation of para-crystalline region as shown in Figure 30. For mild shock in [100] direction, PPTA shows very similar behavior as DFT, see chapter 4.2.3 for more details. 70 Figure 30: Paracrystalline transformation of Kevlar during shock loading in [010] direction. 4.5. Conclusion In conclusion, we used ab-initio molecular dynamics and reactive molecular dynamics to study the shock loading of PPTA based on the multiscale shock technique. The results revealed a highly anisotropic shock response involving a hydrogen bond preserving structural phase transformation, for shock along [100], and polycrystalline Kevlar formation, for shock along the [010] direction. These two distinct shock release mechanisms may help in the understanding of the outstanding performance displayed by aramid fibers and their failure mechanisms, which involve a fibrillation failure process. Our results suggest that the strengthening of aramid polymers produced by the drawing process, which results in highly aligned crystalline domains along the fiber axis, could be significantly enhanced if the superior shock performance of the PPTA [100] direction could be further explored in novel fiber processing and fabric design. 71 5. Multiple Reaction Pathway in Shocked 2,4,6- Triamino-1,3,5-Trinitrobenzene crystal 5.1. Introduction 2,4,6-triamino-1,3,5-trinitrobenzene (TATB, C6H6O6N6), molecular structure shown in Figure 31a and crystal structure shown in Figure 31b, is important HE material due to its technological application.53 It exhibits an unusually high thermal and shock resistances.54 Understanding of its chemical reaction during detonation process with atomistic details may provide an insight about the chemical stability of TATB and eventually leads to the discovery of safer and more powerful HE materials. Thus, both thermal and shock stability of TATB have been studied in great detail by both experiment and theory.48,153–157 In a TATB molecule, the C-NO2 bond dissociation energy is much lower than the C-NH2 bond energy.155 Thus, loss of nitro groups via various possible mechanisms was suggested as primary Figure 31: (a) Single TATB molecule where cyan, grey, red and blue represent carbon, hydrogen, oxygen and nitrogen, respectively. (b) TATB crystalline unit cell. Each unit cell consists of 2 TATB molecules. The unit cell has the lattice parameters of a = 9.010 Å, b = 9.028 Å, c = 6.812 Å, = 108.58˚, = 91.82˚ and = 119.97˚.158 72 pathways. Sharma et al. observed the loss of NO2 via scission of C-NO2 bonds as the primary pathway during thermal and photoinduced decompositions.157 Shock and impact studies by Sharma et al. suggested the formation of benzofurazan or benzofuroxan intermediates.159 First principles electronic structure calculations by Wu et al. suggested that intramolecular hydrogen bonding plays a crucial role in the thermal decomposition of TATB. Hydrogen transfer mediated ring closure leads to the formation of benzofurazan or benzofuroxan intermediates.155 Simultaneous thermogravimetric modulated beam mass spectrometry and time-of-flight velocity spectra measurements also observed the formation of benzofurazan derivatives.54 In contrast to these studies, other spectroscopic and thermal studies by Makashir and Kurian suggested C-NH2 scission as the primary reaction step.160 Farber and Srivastava suggested a carbon-ring cleavage mechanism for the thermal decomposition of TATB.156 Another study using time-of-flight mass spectroscopy by Östmark suggested the formation of nitroso compounds in sub-detonation of TATB.154 On the contrary, reactive force field simulation for the thermal decomposition of TATB suggested the presence of more than one mechanisms.161 Shock decomposition pathways of crystalline TATB were studied by Manna et al. using molecular dynamics (MD) simulations based on the density functional tight binding (SCC-DFTB) method, which also confirmed the presence of nitrogen rich mono and dibenzofurazans heterocycles.153 Heterocyclic clusters formed during the reaction are dynamically stable, and eventually form N2 and CO2 (or CO). Time scales for the formation of these final compounds are in subnanosecond due to the formation of these heterocyclic compounds. These clusters act as reactivity retardants, leading to the shock insensitivity of TATB. However, the formation of these heterocyclic clusters under shock and impact is not well understood, especially within complex heterogeneous shock structures. We study shock-induced decomposition of TATB with atomistic detail. Simulation results provide 73 understanding about the controversial C-NH2 and C-NO2 bond scissions, as well as the formation of N2. We have identified intermolecular pathways for the formation of N2 and H2O. We have also observed the formation of large carbon cluster with high nitrogen and oxygen contents. 5.2. Simulation details The initial system, shown in Figure 32, is composed of 100×4×6 TATB crystalline unit cells in a triclinic MD simulation box of dimensions 139.4 nm × 3.568 nm × 3.870 nm along the a, b and c crystallographic axes. In addition, the MD box includes an empty space of 43 nm at the left end in the a direction in order to observe the evolution of the system over time. In our simulation cell, we have applied a 6.3 nm thick reflective wall at the right end. The total number of TATB molecules is 4,800, amounting to 115,200 atoms. The atomic configuration is first relaxed for 6 ps in the canonical (NVT) ensemble at a temperature of 300 K. After thermalization, we introduce a repulsive flyer at the left end of the system with an initial velocity of 4 km/s for shock loading. The repulsive flyer and the reflective wall at the right end of the system as shown in Figure 32. The interaction between the flyer and TATB molecules is purely repulsive. Figure 32: Schematic diagram of system setup. The orange-colored block at the left end of the system represents a repulsive flyer, while the right end is held by a reflective wall. 74 5.3. Results 5.3.1. Temperature profile At time t = 0, the flyer starts to move towards the TATB crystal and compresses it, while the repulsion between the flyer and TATB decelerates the flyer. After compression for 10 ps, the repulsive flyer starts to retreat. In the compressed TATB crystal, high temperature and pressure trigger exothermic reactions via chemical decomposition. Supersonic detonation wave from decomposition of TATB travels through the crystal with a speed of 9.98±0.5 km/s. Shock wave propagation in a medium can be characterized by a discontinuous change in pressure, temperature or density of the medium. From the temperature profile of the system, we compute the shock velocity for the 1,3,5-triamino-2,4,6-trinitrobenzene (TATB) crystal. Figure 33 shows the time evolution of the temperature discontinuity in the system. Speed of Shock wave can be calculated by temperature front evolution with respected to time. Shock speed in our system is 9.98±0.5 km/s. Figure 33: Temperature profile of a shocked TATB crystal at different times. A discontinuous temperature jump is observed between the unshocked (blue) and shocked (red) regions. 75 Earlier calculation by Manna et al. suggests that hotspot chemistry initiate between 8 and 9 km/s.153 Reaction products jet out towards left as shown in Figure 32. Figure 34 shows the temperature profile as a function of the position along the a-axis and time. As the flyer hits the TATB crystal, the temperature rises rapidly and a shock wave travels through the TATB crystal. The sharp interface between the unreacted (colored blue) and fully/partially reacted (green/red) regions signifies a shock wave passing through the system. The shock wave reaches the end of the system in ~10 ps and the right end of the system starts to heat up due to a rarefaction wave. The fully reacted system expands toward the left empty space as shown by the leftmost red region in figure 3.4. Figure 34: Temperature profile of the system as a function of the position along the a-axis and time. 5.3.2. Fragment Analysis To study chemical reaction pathways, we perform a fragment analysis. A bond is defined by unique cutoff value of bond order for each atom pair.76 A cluster of atoms, connected with 76 covalent bonds, constitutes a molecule. Figure 35 shows the number of molecular products formed as a function of time. Final products include N2, H2O, CO2 and C as shown in Eq. (5.1).162 C 6 H 6 N 6 O 6 → 3H 2 O+3N 2 + 3 2 CO 2 + 9 2 C (5.1) H2O is the first product to form. After that, we observe the rapid production of both N2 and H2O. Unlike in similar HEs such as RDX and HMX, CO2 is the dominant product compared to CO.163 We observe NH3 formation at the initial phase. However, the number of NH3 fragments remains nearly constant after the initial transient. Therefore, NH3 is not included as part of the final products. Similar behavior has been observed during thermal decomposition of TATB by Quenneville et al.162 Figure 35: Number of molecular fragments formed as a function of time for shock speed 9.98 km/s. Red, blue, black and green curves correspond to the formation of N2, H2O, NH3 and CO2, respectively. 77 Time constant for the formation of a product can be estimated by fitting the corresponding number of molecular fragments in Figure 35 to an exponential form as a function of time by using balanced reaction formulae, 𝑛 𝑖 (𝑡 )=𝑛 𝑖∞ [1−exp(− 𝑡 𝜏 )]. (3.2) Here, ni(t) is the number of the ith molecular product at time t, ni∞ is the total number of products formed asymptotically (t → ∞), and τ is the time constant. Least-square fitting provides𝜏 𝐻 2 𝑂 = 750±200 ps,𝜏 𝑁 2 = 500±100 ps and 𝜏 𝐶𝑂 2 = 1,200±500 ps. Similar time scale was suggested by Manna et al.153 These results suggest a single-stage reaction for TATB shock decomposition unlike for RDX. Previous studies by Li et al. for RDX shock decomposition found a two-stage reaction, i.e., rapid production of N2 (𝜏 𝑁 2 = 10±2 ps) and H2O (𝜏 𝐻 2 𝑂 = 30±4 ps) followed by much slower release of CO2(t CO = 900±400 ps).163 Shock-induced TATB decomposition does not exhibit any separation of reaction time in the decomposition stage as observed in RDX. Product formation in TATB is also much slower compared to that in RDX, which is consistent that TATB is very insensitive. This correlation between slow product formation and insensitivity in TATB has been discussed in previous literature, where slow decomposition was attributed to the formation of nitrogen rich heterocyclic compounds. These heterocyclic compounds persist even at very high temperature (3500 K).153 To understand the reaction mechanisms for the formation of the final products, we trace back the origin of constituent atoms of N2 and H2O molecules. Figure 36a and Figure 36b respectively. In Figure 36a, the red line shows the fraction of N2 molecules produced from one TATB molecule, i.e., via intramolecular pathways. The blue line is the fraction of N2 molecules 78 that are formed from more than one TATB molecules, i.e., via intermolecular pathways. We see that N2 molecules are mostly formed via intermolecular pathways. Figure 36: (a) Fraction of N2 molecules that are formed by two distinct pathways as a function of time. The red line shows the fraction of N2 formed by N atoms from a single TATB molecule (intramolecular pathway). The blue line shows N2 originating from more than one TATB molecule (intermolecular pathway). (b) Fraction of H2O molecules that are formed by two distinct pathways as a function of time. The red line shows the fraction of H2O formed by H and O atoms from a single TATB molecule (intramolecular pathway). The blue line shows H2O originating from more than one TATB molecule (intermolecular pathway) In Figure 36b, the red line shows the fraction of H2O products that are formed by H and O atoms that originate from a single TATB molecule, i.e., intramolecular pathways. Also shown by the blue line is the fraction of H2O molecules that are composed of atoms originating from more than two TATB molecules, i.e., intemolecular pathways. We see that H2O products are mostly formed from two adjacent TATB molecules instead of a single TATB molecule. Spatial proximity allows rapid production of H2O. Furthermore, H atoms can react with O atoms located in remote 79 TATB molecules due to their large diffusivity. As a result, intermolecular pathways dominate intramolecular pathway. Similar pathways for the formation of H2O and N2 were suggested by Manna et al.153 Intermolecular pathways seem to be unique for TATB. Similar studies on RDX suggested that 75% of N2 are formed by intramolecular pathways, i.e., N2 is formed from a single RDX molecule. Such a pathway exists in RDX due to the presence of N-NO2 bonds within a RDX molecule, which is absent in a TATB molecule. Figure 37: Intermediate species formed as a function of time. NO2, NO, HONO, OH and NH2 are represented by blue, green, red, black and magenta, respectively. Intermediate species analysis can provide further insight into the reaction mechanisms. Figure 37 shows the number of intermediate fragments as a function of time. We observe the presence of NO2 and NH2 fragments, which can form from the scission of C-NO2 and C-NH2 bonds, respectively. Presence of HONO and NO2 fragments suggests the formation of N2 from C-NO2 bond scission. Similar fragments were also observed during the thermal decomposition of nitromethane.164 This pathway has been suggested for the formation of N2 in RDX and HMX.163,165,166 NH2 from C- NH2 bond scission can form ammonia as final product by abstracting 80 another hydrogen from N-H bond or form N2 by recombination with other NH2 or NO2 group. Experimental studies of high-temperature thermal decomposition of TATB show the formation of NH2. 5.3.3. Detection of new reaction pathway To further understand the mechanism behind the formation of N2, we trace back the initial environment of constituent N atoms. We count the number of N2 molecules formed from different possible combinations, i.e., NH2-NH2, NO2-NH2 and NO2-NO2 pairs. These three pairs represent three different intermolecular reaction pathways I, II and III, respectively. Table 2 shows the percentage of N2 molecules that are formed from different reaction pathways after 40 ps. Surprisingly, we observe that 50% of N2 is formed from NO2-NH2 pairs during shock decomposition. The high percentage of N2 formation from NO2-NH2 pairs may be due to spatial proximity of NH2 and NO2 groups in crystal. In TATB crystal, a NH2 group from a TATB molecule is adjacent to a NO2 group from another TATB molecule as shown in Figure 38 As we compress TATB crystal in the a crystallographic direction, the distance between NH2 and NO2 pairs decreases drastically. This may be the reason for reaction pathway II via NO2-NH2 to be the dominant reaction pathway. Depending on the shock direction in TATB crystal, reactivity and the primary chemical reaction may change.56 In HMX and RDX, N2 products are formed from NO2- NO2 pairs. However, N2 formation from NO2-NO2 pair is very low in the case of TATB. While, it is clear that NO2-NH2 pair dominantly produce N2 from Table 2, the other reaction pathways cannot be ignored quantitatively. Namely, reaction pathways I and III via NH2-NH2 and NO2-NO2 pairs produce the remaining 50% of N2. 81 Figure 38: a 2D layer in the a-b crystallographic plane of TATB crystal. Cyan, grey, red and blue represent carbon, hydrogen, oxygen and nitrogen, respectively. Black circles enclose proximate intermolecular NO2-NH2 pairs In fact, the first step in the shock decomposition of TATB crystal is compressive formation of large clusters from the coalescence of multiple TATB molecules. During the coalescence process, H2O molecules are released. Due to the extreme compression, we observe the distortion of the planar benzene rings in all three pathways. Detailed analyses suggest that reaction pathways diverge to different reaction pathways to form N2 from different fragment pairs. Table 2: Percentage of N2 formed from different fragment pair after 40 ps Fragment pair % of N2 formed NH2-NH2 28 82 NO2-NH2 50 NO2-NO2 22 Figure 39: Time evolution of the fraction of N2 formed from different reaction pathways. Reaction pathway I (from NH2-NH2 pair) to form N2 undergoes transformation via C-NH2 bond scission from two different TATB molecules. After distortion of the planar structure, C-NH2 bond breaks to form two NH2 species. Subsequently, the two NH2 groups are combined to form N2. C-NH2 bond breaking has been observed experimentally under high static pressure. During this process, hydrogen is abstracted by OH or other NH2 species to form H2O or NH3, respectively. Dissociation energy for C-NH2 bond is very high compared to that of C-NO2 bond,155 but spectroscopic studies suggests that hydrogen bonding with C-NO2 group weakens C-NH2 bond.160,167 Presence of NH3 has been experimentally observed,168 which suggests the presence of this pathway. Even though it is not a primary reaction pathway, it accounts for 28% of the production of N2. Reaction pathway II (from NH2-NO2 pair) is found to be the primary channel for the production of N2, accounting for 50% of the total N2 produced. First step in this channel is the loss of oxygen assisted by intramolecular/intermolecular NH2 group. Oxygen loss was also reported experimentally in decomposition of TATB molecules.160 NO species react with nearby 83 NH2 groups from other TATB molecules to form N2 and H2O/OH. Hydrogen and oxygen are abstracted by OH and NH2 species to form water and ammonia. This reaction pathway was also observed in quantum molecular dynamics studies performed by He et al.169Reaction pathway III (from NO2-NO2 pair) contributes only 22% in the formation of N2. This reaction starts with the breaking of C-NO2 bond from two nearby TATB molecules. We also observe the presence of NO2 as an intermediate in shock decomposition of TATB. Two NO2 species form N2 as a final product. C-NO2 bond breaking has been observed experimentally under high static pressure.170 To understand the timeline of these reactions, Figure 39 shows the fraction of N2 formed from different reaction pathways as a function of time. The Figure 39 shows that the NH2-NO2 pathway (reaction pathway II) dominates over other two pathways (reaction pathway I and III). Reaction pathway III dominate over reaction pathway I in the early stage of reaction, since the C-NO2 bond energy is lower than the C-NH2 bond energy. Thus, N2 formation via reaction pathway III is higher. However, after 3 ps, reaction pathway I become the most favorable reaction pathway. Experimental studies of photoinduced decomposition of TATB at ambient and high pressures by Glascoe et al. suggest that elevated pressure may inhibit the hemolysis of C-NO2 bond and block the primary decomposition process.14 Finally, we perform the stoichiometric analysis for the large clusters. Figure 40(a) shows the fraction of the number of large clusters with respect to the total number of atoms. Here, a large cluster is defined as a cluster that contains more than 24 atoms. The number of large clusters increases rapidly first as the flyer starts to compress TATB crystal. When the rarefaction wave returns, the fraction of large clusters increases again. After 20 ps, number of large clusters starts to decline. 84 Figure 40: Time evolution of large clusters (the number of atoms in the cluster, n > 24). (a) Fraction of the number of atoms in large clusters to the total number of atoms. (b) Atomic composition of the large clusters as a function of time. Red and blue lines Figure 40(b) shows the carbon-to-oxygen (red line) and carbon-to-nitrogen (blue line) ratios in the large clusters. We observe that large clusters are nitrogen and oxygen rich, which delay the formation of stable products such as N2, CO2 and H2O in later stages. Most of the final products seep out of these large clusters. Quantum molecular dynamics studies by Manna et al. found the presence of nitrogen rich clusters, which impede the formation of N2.153 Oxygen rich cluster were also reported in shock simulation of RDX163 and thermal decomposition of TATB.161 5.4. Conclusion In summary, our RMD simulations reveal that N2 production in shock decomposition of TATB is a multi-channel process. Reaction mechanism studies suggest that a reaction pathway via NH2-NO2 pairs is the primary reaction pathway for the formation of N2. However, other two pathways (from NH2-NH2 and NO2-NO2 pairs) also participate significantly in the formation of N2. Primary reaction stages include the formation of N2 and H2O rapidly and then somewhat slower formation of CO2. Most of the products formed via intermolecular pathways 85 6. Electric field induced spontaneous polarization switching in Al(1-x)ScxN 6.1. Introduction Spontaneous polarization switching in ferroelectric has been leveraged into RF devices, non- volatile memory, integrated circuit technology, sensors, actuators and neuromorphic computing.58– 60 Traditional ferroelectrics, such as HfO2/ZrO2 and perovskites, presented many technological challenges due to their long and diffusive transition in hysteresis curve between two polarization states.59,61 Recently discovered Nitride ferroelectrics, Al1-xScxN (x = 0.27 – 0.43), counter such issue and show a sharp square hysteresis curves even in polycrystalline films.62 Nitride ferroelectrics, Al1-xScxN (x = 0.27 – 0.43), are highly attractive for device purposes as their spontaneous polarization large (60—100 µC/cm²), tunable coercive voltage (180-400 V/𝜇 m), and low fabrication temperature(300K) compared to tradition ferroelectric such as PZT and Fe- HfO2.62,63 Nitride ferroelectric have opened are experimentally produced as polycrystalline bulk materials.62,171 The local structure of Al1-xScxN and role of Sc in tunability of ferroelectrice is still unknown. Further, the process leading to extremely sharp transition between two state in states are also unknown. In Chapter 6.2, we will perform quantum molecular dynamics to create Al1-xScxN alloys with different concentration and study the local structure of these alloys. In chapter 6.3, using a Berry-phase approach64–66 with QMD we will model the electric field induced polarization switching in Nitride ferroelectrics. Chapter 6.4 shows the role of Sc concentration in Al1-xScxN alloys. 86 6.2. Computational synthesis of Al(1-x)ScxN alloys 6.2.1. Local structure of Al(1-x)ScxN alloys Al1-xScxN alloys are created by adding AlN and ScN in different concentration. Mixture of both crystals were heated to 4000K temperature for 10 ps in the NPT ensemble. The system is subsequently cooled to 300 K temperature in 10 ps. The cycle of melt quench was performed for 5 times. After melt quench, we perform optimization of the system with using BFGS method while relaxing both lattice constants and atomic positions. Figure 41a and b the local structure of Al0.67Sc0.33N and Al0.60Sc0.40N. We observe the formation of 4 coordinated Sc atoms in lower concentration of ScN suggesting the presence of ScN in wurtzite structure, shown in Figure 41a. Potential energy surface calculation shows that ScN in wurtzite structure is extremely unstable172. However, in presence of low concentration of Sc, we observe the formation of wurtzite structure. Figure 41: local structure of Al1-xScxN ferroelectrics. a) Al0.67Sc0.33N local structure b) Al0.60Sc0.40N local structure. Cyan, black and brown spheres correspond to nitrogen, aluminum and scandium atoms. Red and blue sphere highlight ScN local structure. 87 Figure 41b shows the local structure of Al0.60Sc0.40N. Increasing the concentration of Sc leads to the formation of six coordinated Sc atoms, as observed in ScN rocksalt structure. The local structure suggests that introduction of Sc in AlN creates a metastable state of Al1-xScxN alloys with point group symmetry of wurtzite structure. The AlN ground state is destabilized by Sc leading to a lowered switching barrier of between two minima of piezoelectric and thus converting it into ferroelectric. 6.2.2. Phase diagram of Al(1-x)ScxN alloys Figure 42: Phase diagram of Al(1-x)ScxN alloys Figure 42 shows the phase diagram of Al(1-x)ScxN alloys. The red and blue curve corresponds to the rocksalt and wurtzite structure of alloys, respectively. We observe a transition point at 0.43 corresponding to the switching two aforementioned phases. Below 0.43, we observe wurtzite structure is dominant while above 0.43, rock-salt structure is dominant. The ferroelectricity in nitride alloys is observed experimentally between 𝑥 =0.27-0.43, highlighted by 88 the yellow region.62 The coercive field required for transition between two minima to becomes higher for concentration below 0.27 leading to formation of piezoelectric alloys. On the contrary, rocksalt structure is stabilized above 0.43 with no spontaneous polarization. The alloys acts as dielectric materials. 6.3. Switching pathway of Al(1-x)ScxN alloys The models used for switching simulations along the crystallographic directions [001], we created systems with 256 atoms. We have created two system corresponding to Al0.75Sc0.25N and Al0.5Sc0.5N. Initially, the system is relaxed using a Quasi-Newton method (BFGS algorithm) at 0 GPa pressure. After relaxation, the finite electric field is applied with berry phase approach outlined by Umari et al within periodic boundary condition.64 Time-step for all the simulation is 50 a.u. 6.3.1. Energy barrier calculation Figure 43a shown the crystal structure of Al0.75Sc0.25N initial configuration. Here nitrogen atoms are in down position (down and up position are arbitrarily defined). Figure 43b shows the configuration of switched structure in presence of electric field. The switching is performed between two state at temperature T= 300K in presence of electric field 3500 V/𝜇𝑚 . It is important to note that the switching field is 8 times higher than the experimental value.62 Al0.5Sc0.5N alloys also undergoes switching via same pathway (figure not shown). Figure 43c shows the potential energy surface for switching pathway for both Al0.75Sc0.25N (red curve) and Al0.5Sc0.5N (blue curve). Switching barrier for Al0.75Sc0.25N and Al0.5Sc0.5N is 0.29 eV and 0.25 eV. From the MD trajectories we have computed the average charges and displacement in N atoms for Al0.75Sc0.25N which is -0.85 and -0.96. Using the aforementioned switching value, electric field required for switching between two polarization is 355 V/𝜇𝑚 , which is in good agreement with the experiment 89 (coercive field= 380 V/𝜇𝑚 .62 Similarly switching field for Al0.5Sc0.5N alloys is 264 V/𝜇𝑚 . Further from Figure 43c, it is clear that addition of Sc lowers the coercive field required the spontaneous polarization switching. Figure 43: Al0.75Sc0.25N configuration and energy barrier. Al0.75Sc0.25N in down (a)/up (b) position. (c) energy barrier for switching of Al0.75Sc0.25N (red curve) and Al0.75Sc0.25N (blue curve). Grey and yellow spheres correspond to Al and Sc atoms respectively. Blue/red sphere corresponds to Ndown/Nup. 6.3.2. Switching Metastable Layered nitride Previously, it has been hypothesized that the switching in Al(1-x)ScxN alloys occurs via metastable layered nitride.173 To further understand the switching pathway of spontaneous polarization, we plot position of Al (red curves) and Sc (blue curves) as a function of time for one layer in Figure 44. From the figure, we do not observe a bulk switching in Al0.75Sc0.25N alloys. Each Al atoms switch at different time suggesting a nucleation of switching and its growth in the system. Most of the Sc atoms follow the same trajectory. Further, we observe an abrupt change in positions of Al (~ 0.03 ps), while change in Sc positions occurs 0.3 ps. This suggests that Al does not goes through metastable 5 coordinated phase. However, we observe the formation of such 5 90 coordinated ScN is observed for Sc atoms. Snapshot in Figure 44 shows the 5 coordinated Sc formation. Figure 44: position of Al (red curves) and Sc (blue curves) as a function of time for one layer. Snapshot of Sc and its surrounding N atoms at different time are shown by around image. Yellow and blue sphere corresponds to the Sc and N, respectively. Formation of h-ScN (5 coordinated Sc) in Al(1-x)ScxN alloys has wide implication. h-ScN is a metastable phase compared to the wurtzite structure.172 Thus, increasing concentration of Sc in wurtzite structure will make the system increasingly unstable while stabilizing the transition point during switching. This allows a lower switching field as observed both experimentally62 and theoretically in this work. Second, the homogeneous distribution of Sc in Al(1-x)ScxN will show higher barrier compared to heterogenous distribution due to ease of formation of h-ScN in presence of higher local concentration of Sc. Thus, a melt-quench process should provide a better ferroelectric. 91 6.4. Conclusion In conclusion, here we have used adiabatic quantum molecular dynamics to perform a computational synthesis of Al(1-x)ScxN. We modeled the polarization switching of nitride ferroelectric induced by electric field. Our simulation results are in excellent agreement with experiment. We have also computed the switching barrier and switching pathway. Here, we propose a mixed pathway contrary to previously hypothesized, where only Sc atoms form nitride hexagonal phase. Our study should allow us to perform high throughput computational screening of ferroelectrics. 92 7. Machine Learning Assisted Principle of Least Action Formalism for Schrodinger Equation 7.1. Introduction The principle of least action (PLA) is the cornerstone of physics. It forms the basis for Newtonian mechanics, special and general theory of relativity, and electrodynamics. Fermat developed a slight variant of PLA for wave optics and Feynman formulated the PLA for quantum mechanics in his Ph.D. thesis.174 The basic idea of PLA is to find the path that minimizes the action integral. In mechanics, the action integral A is the difference between time-averaged kinetic and potential energies between the initial and final configurations at two different instants of time: 𝐴 = ∫[𝐾 ( 𝑞 ̇(𝑡 ))−𝑉 (𝑞 (𝑡 ))]𝑑𝑡 (7.1) where the integrand is the Lagrangian of the system, the potential energy V is a function of generalized coordinates q(t), and the kinetic energy K depends on generalized velocities 𝑞 ̇(𝑡 ) . Minimization of the action integral gives Euler-Lagrange equations, which form the basis of molecular dynamics (MD) simulations. The PLA can be formulated for various kinds of physical problems. For example, the Poisson’s equation, ∇ 2 𝜙 = −𝜌 /𝜖 0 (7.2) can be recast as, 𝑈 𝑃𝐸 ∗ = 𝜖 0 2 ∫ (∇𝜙 ) 2 𝑑 Ω− ∫ 𝜌𝜙𝑑 Ω (7.3) where 𝜙 is the electrostatic potential, 𝜌 is the charge density, 𝜖 0 is the dielectric constant, and is the volume of the system. The Poisson equation follows from the minimization of 𝑈 𝑃𝐸 ∗ with respect to 𝜙 .175 Similarly, the time-independent Schrödinger equation for a quantum particle 93 (− ℏ 2 2𝑚 ∇ 2 +𝑉 (𝒓 ))𝜓 (𝒓 )=𝐸𝜓 (𝒓 ) (7.4) can be recast as 𝑈 𝑆𝐸 ∗ = ℏ 2 2𝑚 ∫(∇𝜓 ) 2 𝑑 3 𝒓 +∫ 𝜓 ∗ (𝑉 (𝒓 )−𝐸 )𝜓 𝑑 3 𝒓 (7.5) where 𝑚 is the mass, 𝑉 (𝒓 ) is the potential energy, and 𝐸 is the total energy of the particle. The Schrödinger equation can be obtained by minimizing 𝑈 𝑆𝐸 ∗ with respect to the wavefunction 𝜓 . In this paper, we examine whether a neural network can learn to solve Schrödinger equation by minimizing 𝑈 𝑆𝐸 ∗ . Further, we use machine learning (ML) to solve the Schrödinger equation directly with the automatic differentiation (autodiff) method.71 In recent years, significant progress has been made in solving ordinary differential equations (ODEs) and partial differential equations (PDEs) with autodiff and physics informed neural networks. Raissi et al. have attempted to solve nonlinear PDEs of the form,176 𝑢 𝑡 +𝒩 [𝑢 ;𝜆 ]=0, 𝑥 ∈Ω,𝑡 ∈[0.𝑇 ] by employing deep neural networks as universal function approximators. Here 𝒩 (𝑢 ;𝜆 ) is a nonlinear operator parameterized by 𝜆 , is a subset of ℝ 𝐷 , and 𝑢 (𝑡 ,𝑥 ) is the latent deep learning solution for the PDE. Raissi et al. have used both continuous and discreet time models to solve various non-linear problems such as the time-dependent Schrödinger equation, Burger’s equation and Navier-Strokes equation. They solve the time-dependent Schrödinger equation with a continuous time model using deep learning with a loss function derived from the physical equation itself. It is important to note that parameters are shared between 𝑢 (𝑡 ,𝑥 ) and 𝒩 (𝑢 (𝑡 );𝜆 ) and the loss function includes constraints arising from initial and final boundary conditions. Their loss functions consist of mean square errors, 𝑀𝑆𝐸 =𝑀𝑆 𝐸 𝑢 +𝑀𝑆 𝐸 𝑓 (7.7) 94 where 𝑀𝑆 𝐸 𝑢 = 1 𝑁 𝑢 ∑|𝑢 (𝑡 𝑢 𝑖 ,𝑥 𝑢 𝑖 )−𝑢 𝑖 | 2 𝑁 𝑢 𝑖=1 (7.7a) and 𝑀𝑆 𝐸 𝑓 = 1 𝑁 𝑢 ∑ |𝑓 (𝑡 𝑓 𝑖 ,𝑥 𝑓 𝑖 )| 2 𝑁 𝑢 𝑖=1 ; 𝑓 ≔𝑢 𝑡 +𝑁 [𝑢 ] (7.7b) In this paper, we show that a feed-forward neural network is capable of learning quantum mechanics via the PLA. We have chosen hydrogen atom as a model for the neural network. We find that the ground and excited state wavefunctions calculated by the neural network are in close agreement with the exact analytical results. We have also used autodiff in conjunction with a self- consistent scheme to solve the Schrodinger equation for the hydrogen atom. This approach provides energies and wavefunctions of ground and excited states that are in excellent agreement with the analytical results. Remarkably, the combined ML and self-consistent methods provide exact wave functions with the Gram-Schmidt orthogonalization. 7.2. Method In this work, we have used a feed-forward neural network to minimize the principle of least action for the radial part of a hydrogen atom subject to the constraint arising from the boundary conditions and wave-function orthogonalization. The radial Schrödinger equation (S.E.) for the hydrogen atom is (− 1 2𝑟 𝜕 2 𝜕 𝑟 2 + 𝑙(𝑙 +1) 2𝑟 2 − 1 𝑟 𝜓 )=𝐸𝜓 (7.8) which follows from the minimization of, 95 𝑈 𝐻 = 1 2 ∫ ( 𝜕𝜓 𝜕𝑟 ) 2 𝑟 0 𝑑𝑟 +∫ 𝜓 ∗ 𝑟 0 ( 𝑙(𝑙 +1) 2𝑟 2 − 1 𝑟 −𝐸 )𝜓𝑑𝑟 (7.9) with respect to the wave function. Here 𝑟 is the distance between the nucleus and the electron, and 𝜓 and E are the eigen state and the corresponding eigenvalue, respectively. The architecture of the neural network we have used to minimize 𝑈 𝐻 is shown in Figure 45. The input layer stores the spatial grid on which the wavefunction is represented. Each grid point in the input layer is connected to every unit in the hidden layer. Typically, the hidden layer has 20 or 40 units. The data in the input layer undergoes linear transformation Wr + b, where W is the weight matrix connecting input-layer and hidden-layer units and b is a vector representing the hidden-unit biases. The linear transformation is followed by a non-linear transformation with the tanh activation function. The output of the hidden layer is the radial wavefunction. The loss function is the total mean square error (MSE) consisting of 𝑀𝑆𝐸 = 𝑈 𝐻 2 +𝑀𝑆 𝐸 𝑢 (7.10) where 𝑀𝑆 𝐸 𝑢 contains constraints due to boundary conditions and wave-function orthogonalization: 𝑀𝑆 𝐸 𝑢 =𝑀𝑆 𝐸 𝑛𝑜𝑟𝑚 +𝑀𝑆 𝐸 𝑟 0 +𝑀𝑆 𝐸 𝑟 ∞ +𝑀𝑆 𝐸 𝑑 ∞ +𝑀𝑆 𝐸 𝑜𝑟𝑡 ℎ𝑜 where 𝑀𝑆 𝐸 𝑛𝑜𝑟𝑚 =(⟨𝜓 |𝜓 ⟩−1) 2 𝑀𝑆 𝐸 𝑟 0 =𝜓 | 𝑟 =0 2 𝑀 𝑆𝐸 𝑟 ∞ =𝜓 | 𝑟 =∞ 2 𝑀𝑆 𝐷 𝑑 ∞ = 𝑑𝜓 𝑑𝑟 | 𝑟 =∞ 2 𝑀𝑆 𝐷 𝑜𝑟𝑡 ℎ𝑜 =∑(⟨𝜓 𝑛𝑠 |𝜓 𝑖𝑠 ⟩−1) 2 𝑛 −1 𝑖 =0 96 The mean square error 𝑀𝑆 𝐸 𝑢 includes constraints due to wave-function orthonormalization (𝑀𝑆 𝐸 𝑛𝑜𝑟𝑚 ) and boundary conditions. The latter, 𝑀𝑆 𝐸 𝑟 0 and 𝑀 𝑆𝐸 𝑟 ∞ , ensure that the wavefunctions vanish at the origin and far away from the origin. In addition, the first derivative of the wave function must also vanish at large distances from the nucleus. Our cutoff for 1s, 2s and 2p orbitals is 30 Å and we use 30,000 grid points equally spaced along the radial direction. The mean square error between the neural network wavefunction and the analytical result becomes negligible after 10000 epochs. Figure 45: A scheme of feed forward neural network employed in our calculation. We have used tanh as a non-linear unit. We have also used machine learning to solve the S.E. for a hydrogen atom. In this case, we minimize the loss function 𝑀𝑆𝐸 = 𝑀𝑆 𝐸 𝑒𝑞 +𝑀𝑆 𝐸 𝑢 (7.11) where 𝑀𝑆 𝐸 𝑒𝑞 =(− 𝑟 2 𝜕 2 𝜓 𝜕 𝑟 2 + 𝑙(𝑙 +1) 2𝑟 𝜓 −(1+𝐸𝑟 )𝜓 ) 2 and 𝑀𝑆 𝐸 𝑢 =𝑀𝑆 𝐸 𝑛𝑜𝑟𝑚 +𝑀𝑆 𝐸 𝑟 0 +𝑀𝑆 𝐸 𝑟 ∞ +𝑀𝑆 𝐸 𝑑 ∞ The loss function for the S.E is modified by multiplying with 𝑟 to avoid error explosion for small values of 𝑟 . Note that here the loss 𝑀𝑆 𝐸 𝑢 does not include the constraint due to wave-function 97 orthogonalization. All other terms in 𝑀𝑆 𝐸 𝑢 are the same as in Eq. 6.10. We use the same cutoff, rc = 30 Å for 1s, 2s and 2p orbitals and rc = 35 Å for 3d orbitals. We use 200 grid points equally spaced along the radial direction. The accuracy of our method is assessed by comparing the radial wavefunction against the analytical solution. We employ autodiff to compute the term 𝑑 2 𝜓 𝑑 𝑟 2 in the loss function. Automatic differentiation is now frequently employed to do backpropagation in deep learning. We use a straightforward self-consistent field (SCF) algorithm to find the stationary eigenstates and the corresponding eigenvalues. The algorithm is given in Algorithm 2. Function SCF(): ========================= For any E: 𝐸 𝑝𝑟𝑒𝑣 =0 ========================= While (𝐸 −𝐸 𝑝𝑟𝑒𝑣 ) <𝛿 : For epoch in epochs Compute MSE [𝜓 ;𝐸 ] Update weight of neural network Compute updated wave function 𝜓 (𝑟 ) Normalized 𝜓 = 𝜓 ⟨𝜓 |𝜓 ⟩ Compute new energy using equation (− 1 2 𝜕 2 𝜓 𝜕 𝑟 2 + 𝑙 (𝑙+1) 2𝑟 2 𝜓 −( 1 𝑟 )𝜓 ) if (𝐸 −𝐸 𝑝𝑟𝑒𝑣 )<𝛿 : return 𝐸 ,𝜓 else: 𝐸 𝑝𝑟𝑒𝑣 =𝐸 continue Algorithm 2: A naïve SCF scheme to compute a stationary eigenstate and its corresponding eigen-energy for hydrogen atoms. In the SCF implementation, the neural network computes the best 𝜓 for an arbitrary value of E. The normalized 𝜓 from the output of the neural network is used to compute the corresponding 98 energy from the S. E. We perform a SCF check over energy and if the check is successful, we return the final eigen value and eigen vector. However, if the SCF check fails, the computed energy 𝐸 is used as an input in the next iteration and the iterations are continued until the energy converges. The final output, 𝜓 and E, after the convergence of SCF iterations are the eigen vector and eigen value of the S.E. for the hydrogen atom. 7.3. Results Figure 46 shows the loss, given by equation 7.10, as a function of epoch for the PLA method. We can observe the after 6000 steps, the error is almost zero. Figure 46b, c and d show 1s, 2s, and 2p radial wave function, respectively, as a function of epoch. It is important to note that even though error is extremely low after 4000 epoch, wavefunction shape do not match with analytical solution. After 10000 epochs, we can observe that analytical, shown by blue curve, and neural network solution, shown by red curve, are in good agreement. Figure 46: Training results of neural network to obtain radial wave function using PLA method. a) loss as a function of epochs b-d) 1s, 2s, and 2p 99 It is important to note that PLA method optimizes the action integral to find the optimal wavefunction. However, it fails to find excited state solution without Gram-Schmidt orthogonalization method. The orthonormality condition for radial part of wavefunction for hydrogen atom is defined by ∫ 𝜓 𝑛𝑙 ∗ 𝑟 0 𝜓 𝑛 ′ 𝑙 𝑑 3 𝑟 =𝛿 𝑛 𝑛 ′ (6.9) Where 𝑛 ,𝑛 ′ are principle quantum numbers and 𝑙 is azimuthal quantum number. Figure 47: Training results of neural network to obtain radial wave function using Hamiltonian method. a) loss as a function of epochs b-d) 1s, 2s, and 3s Figure 46 shows the loss, given by equation 7.10, as a function of epoch for the PLA method. We can observe that after 6000 steps, the error is almost zero. Figure 46b, c and d show 1s, 2s, and 2p radial wave function, respectively, as a function of epoch. It is important to note that even though error is extremely low after 4000 epoch, wavefunction shape do not match with analytical 100 solution. After 10000 epochs, we can observe that analytical, shown by blue curve, and neural network solution, shown by red curve, are in good agreement. Figure 47 shows the loss, given by equation 7.11, as a function of epoch for the Hamiltonian method. The number of epochs required for this method far exceed the number of epochs required for PLA method. Figure 47b, c and d show 1s, 2s, and 3s radial wave function, respectively, as a function of epoch. We can observe that as the number of nodes in the wavefunction increase, the neural network fitting becomes worse due to complexity of the curve. However, the final wave function and energy shows excellent agreement with the literature. Solution of Hamiltonian method does not require the Gram-Schmidt orthogonalization method. Figure 48: SCF iteration and resulting 1s, 2s, 2p orbitals So far, we have discussed the wavefunction solution where energy is parameter for loss function. However, we use a straightforward self-consistent field (SCF) algorithm to find the stationary eigenstates and the corresponding eigenvalues. Figure 48 show the result of SCF algorithm where we provide arbitrary value of energy into Algorithm 2. The energy and wavefunction converge to 101 nearest stationary state. It is important to not again that excited states are determined without the application of Gram-Schmidt orthogonalization method 7.4. Conclusion Here, we have also used autodiff in conjunction with a self-consistent scheme to solve the Schrodinger equation for the hydrogen atom. This approach provides energies and wavefunctions of ground and excited states that are in excellent agreement with the analytical results. Remarkably, the combined ML and self-consistent methods provide exact wave functions with the Gram- Schmidt orthogonalization. 102 8. Conclusion This thesis has explored modeling of quantum materials under non-equilibrium condition using quantum molecular dynamics—specifically electric field, photoexcitation and shock. These conditions are prevalent in the experimental conditions and are often used to study the quantum material’s exotic properties. Modeling their dynamic response with accurate spatial and temporal resolution can be achieved by non-equilibrium quantum molecular dynamics simulations. First, we studied the photoexcitation response of GeTe and Sb2Te3, a common phase change material, to find the non-thermal amorphization pathway. Phase change materials are of high interest for low-power, high-throughput storage devices for next generation neuromorphic computing technologies. Our results highlight an energy-efficient and ultrafast amorphization pathway that could be used to enhance the performance of phase change materials based opto- electronic devices. Second, we modeled shock in aramid fibers with quantum and classical molecular dynamics simulations using multiscale shock theory. Aramid fibers are widely used as lightweight, shock-resistance, reinforcing materials. We found an interplay between hydrogen bond preserving and non-preserving shock release mechanisms are critical to understanding the shock performance of aramid fibers. Using large-scale classical molecular dynamics simulations, we further found the formation of paracrystalline Kevlar. This process was impossible to study by quantum molecular dynamics due to high computational cost, suggesting the importance of multiscale simulation approach to study these non-equilibrium materials processes encompassing large spatiotemporal scales. 103 Third, we modeled non-equilibrium shock simulation to understand the detonation process in 2,4,6-triamino-1,3,5-trinitrobenzene (TATB), where we revealed the existence of three competing intermolecular pathways for the formation of N2. Finally, we used adiabatic quantum molecular dynamics combined with a Berry-phase approach to perform computational synthesis of Nitride ferroelectrics, Al1-xScxN (x = 0.27 – 0.43), thereby modeling electric field-induced switching of polarization. Our simulation results are in excellent agreement with experiment, while identifying the switching pathway and significant effect of Sc doping in Al1-xScxN. Apart from all the traditional molecular dynamics techniques, we also used machine learning to extract ground and excited states from Schrödinger equation. machine learning technique used to learn Schrödinger equation provides energies and wavefunctions of ground and excited states that are in excellent agreement with the analytical results. With the coming of post-Moore era, quantum materials have become exceedingly important to explore new paradigm of computing. Quantum molecular dynamics and machine learning can be utilized synergistically to explore the high-dimensional materials space and model fundamental and exotic properties occurring in the materials with high fidelity. In the present era of quantum materials where non-equilibrium processes are prevalent, this thesis has demonstrated a promising roadmap to study non-equilibrium processes using quantum molecular dynamics and machine learning. 104 References (1) Aaronson, S. The Limits of Quantum. Sci. Am. 2007, No. March, 62–69. (2) Arute, F.; Arya, K.; Babbush, R.; Bacon, D.; Bardin, J. C.; Barends, R.; Biswas, R.; Boixo, S.; Brandao, F. G. S. L.; Buell, D. A.; et al. Quantum Supremacy Using a Programmable Superconducting Processor. Nature 2019, 574, 505–510. (3) Calimera, A.; Macii, E.; Poncino, M. The Human Brain Project and Neuromorphic Computing. Funct. Neurol. 2013, 28, 191–196. (4) Samarth, N. Quantum Materials Discovery from a Synthesis Perspective. Nat. Mater. 2017, 16, 1068–1076. (5) Ball, P. Quantum Materials: Where Many Paths Meet. MRS Bull. 2017, 42, 698–705. (6) Buzzi, M.; Först, M.; Mankowsky, R.; Cavalleri, A. Probing Dynamics in Quantum Materials with Femtosecond X-Rays. Nat. Rev. Mater. 2018, 3, 299–311. (7) Graves, C. E.; Reid, A. H.; Wang, T.; Wu, B.; de Jong, S.; Vahaplar, K.; Radu, I.; Bernstein, D. P.; Messerschmidt, M.; Müller, L.; et al. Nanoscale Spin Reversal by Non- Local Angular Momentum Transfer Following Ultrafast Laser Excitation in Ferrimagnetic GdFeCo. Nat. Mater. 2013, 12, 293–298. (8) Stanciu, C. D.; Hansteen, F.; Kimel, A. V.; Kirilyuk, A.; Tsukamoto, A.; Itoh, A.; Rasing, T. All-Optical Magnetic Recording with Circularly Polarized Light. Phys. Rev. Lett. 2007, 99, 047601. (9) Chen, F.; Zhu, Y.; Liu, S.; Qi, Y.; Hwang, H. Y.; Brandt, N. C.; Lu, J.; Quirin, F.; Enquist, H.; Zalden, P.; et al. Ultrafast Terahertz-Field-Driven Ionic Response in Ferroelectric BaTiO3. Phys. Rev. B 2016, 94, 180104. 105 (10) de Jong, S.; Kukreja, R.; Trabant, C.; Pontius, N.; Chang, C. F.; Kachel, T.; Beye, M.; Sorgenfrei, F.; Back, C. H.; Bräuer, B.; et al. Speed Limit of the Insulator–Metal Transition in Magnetite. Nat. Mater. 2013, 12, 882–886. (11) Berruto, G.; Madan, I.; Murooka, Y.; Vanacore, G. M.; Pomarico, E.; Rajeswari, J.; Lamb, R.; Huang, P.; Kruchkov, A. J.; Togawa, Y.; et al. Laser-Induced Skyrmion Writing and Erasing in an Ultrafast Cryo-Lorentz Transmission Electron Microscope. Phys. Rev. Lett. 2018, 120, 117201. (12) Basov, D. N.; Averitt, R. D.; Hsieh, D. Towards Properties on Demand in Quantum Materials. Nat. Mater. 2017, 16, 1077–1088. (13) Hearmon, A. J.; Fabrizi, F.; Chapon, L. C.; Johnson, R. D.; Prabhakaran, D.; Streltsov, S. V.; Brown, P. J.; Radaelli, P. G. Electric Field Control of the Magnetic Chiralities in Ferroaxial Multiferroic RbFe(MoO4)2. Phys. Rev. Lett. 2012, 108, 237201. (14) Glascoe, E. A.; Zaug, J. M.; Armstrong, M. R.; Crowhurst, J. C.; Grant, C. D.; Fried, L. E. Nanosecond Time-Resolved and Steady-State Infrared Studies of Photoinduced Decomposition of TATB at Ambient and Elevated Pressure. J. Phys. Chem. A 2009, 113, 5881–5887. (15) Wu, C.-C.; Jariwala, D.; Sangwan, V. K.; Marks, T. J.; Hersam, M. C.; Lauhon, L. J. Elucidating the Photoresponse of Ultrathin MoS 2 Field-Effect Transistors by Scanning Photocurrent Microscopy. J. Phys. Chem. Lett. 2013, 4, 2508–2513. (16) Wuttig, M.; Yamada, N. Phase-Change Materials for Rewriteable Data Storage. Nat. Mater. 2007, 6, 824–832. (17) Peng, C.; Wu, L.; Song, Z.; Rao, F.; Zhu, M.; Li, X.; Liu, B.; Cheng, L.; Feng, S.; Yang, P.; et al. Performance Improvement of Sb2Te3 Phase Change Material by Al Doping. Appl. 106 Surf. Sci. 2011, 257, 10667–10670. (18) Zier, T.; Zijlstra, E. S.; Kalitsov, A.; Theodonis, I.; Garcia, M. E. Signatures of Nonthermal Melting. Struct. Dyn. 2015, 2, 054101. (19) Ielmini, D.; Lacaita, A. L. Phase Change Materials in Non-Volatile Storage. Mater. Today 2011, 14, 600–607. (20) Hamann, H. F.; O’Boyle, M.; Martin, Y. C.; Rooks, M.; Wickramasinghe, H. K. Ultra- High-Density Phase-Change Storage and Memory. Nat. Mater. 2006, 5, 383–387. (21) Lee, T. H.; Loke, D.; Huang, K.-J.; Wang, W.-J.; Elliott, S. R. Tailoring Transient- Amorphous States: Towards Fast and Power-Efficient Phase-Change Memory and Neuromorphic Computing. Adv. Mater. 2014, 26, 7493–7498. (22) Skelton, J. M.; Loke, D.; Lee, T.; Elliott, S. R. Ab Initio Molecular-Dynamics Simulation of Neuromorphic Computing in Phase-Change Memory Materials. ACS Appl. Mater. Interfaces 2015, 7, 14223–14230. (23) Chatzi, E. G.; Koenig, J. L. Morphology and Structure of Kevlar Fibers: A Review. Polym. Plast. Technol. Eng. 1987, 26, 229–270. (24) Northolt, M. G.; Den Decker, P.; Picken, S. J.; Baltussen, J. J. M. M.; Schlatmann, R. The Tensile Strength of Polymer Fibres. In Advances in Polymer Science; Springer-Verlag: Berlin/Heidelberg, 2005; Vol. 178, pp 1–108. (25) Guoqi, Z.; Goldsmith, W.; Dharan, C. H. Penetration of Laminated Kevlar by Projectiles—I. Experimental Investigation. Int. J. Solids Struct. 1992, 29, 399–420. (26) Mayo, J. B.; Wetzel, E. D. Cut Resistance and Failure of High-Performance Single Fibers. Text. Res. J. 2014, 84, 1233–1246. (27) Northolt, M. G. X-Ray Diffraction Study of Poly(p-Phenylene Terephthalamide) Fibres. 107 Eur. Polym. J. 1974, 10, 799–804. (28) Dobb, M. G.; Hindeleh, A. M.; Johnson, D. J.; Saville, B. P. Lattice Resolution in an Electron-Beam Sensitive Polymer. Nature 1975, 253, 189–190. (29) Young, R. J.; Lu, D.; Day, R. J.; Knoff, W. F.; Davis, H. A. Relationship between Structure and Mechanical Properties for Aramid Fibres. J. Mater. Sci. 1992, 27, 5431– 5440. (30) Liu, J.; Geil, P. H. Crystal Structure and Morphology of Poly(Ethylene Terephthalate) Single Crystals Prepared by Melt Polymerization. J. Macromol. Sci. Part B 1997, 36, 61– 85. (31) Liu, J.; Cheng, S. Z. D.; Geil, P. H. Morphology and Crystal Structure in Single Crystals of Poly(p-Phenylene Terephthalamide) Prepared by Melt Polymerization. Polymer (Guildf). 1996, 37, 1413–1430. (32) Rao, Y.; Waddon, A. J.; Farris, R. J. Structure–Property Relation in Poly(p-Phenylene Terephthalamide) (PPTA) Fibers. Polymer (Guildf). 2001, 42, 5937–5946. (33) Lee, K.-G.; Barton, R.; Schultz, J. M. Structure and Property Development in Poly(p- Phenylene Terephthalamide) during Heat Treatment under Tension. J. Polym. Sci. Part B Polym. Phys. 1995, 33, 1–14. (34) Ahmed, D.; Hongpeng, Z.; Haijuan, K.; Jing, L.; Yu, M.; Muhuo, Y. Microstructural Developments of Poly ( p-Phenylene Terephthalamide ) Fibers During Heat Treatment Process : A Review. Mater. Researh 2014, 17, 1180–1200. (35) O’Connor, T. C.; Robbins, M. O. Chain Ends and the Ultimate Strength of Polyethylene Fibers. ACS Macro Lett. 2016, 5, 263–267. (36) Chowdhury, S.; Sockalingam, S.; Gillespie, J. Molecular Dynamics Modeling of the 108 Effect of Axial and Transverse Compression on the Residual Tensile Properties of Ballistic Fiber. Fibers 2017, 5, 7. (37) Mercer, B.; Zywicz, E.; Papadopoulos, P. Molecular Dynamics Modeling of PPTA Crystallite Mechanical Properties in the Presence of Defects. Polymer (Guildf). 2017, 114, 329–347. (38) Kunisada, T.; Ushiyama, H.; Yamashita, K. A New Implementation of Ab Initio Ehrenfest Dynamics Using Electronic Configuration Basis: Exact Formulation with Molecular Orbital Connection and Effective Propagation Scheme with Locally Quasi-Diabatic Representation. Int. J. Quantum Chem. 2016, 116, 1205–1213. (39) Grujicic, M.; Bell, W. C.; Glomski, P. S.; Pandurangan, B.; Yen, C.-F. F.; Cheeseman, B. A. Filament-Level Modeling of Aramid-Based High-Performance Structural Materials. J. Mater. Eng. Perform. 2011, 20, 1401–1413. (40) Ortiz, V.; Nielsen, S. O.; Discher, D. E.; Klein, M. L.; Lipowsky, R.; Shillcock, J. Dissipative Particle Dynamics Simulations of Polymersomes. J. Phys. Chem. B 2005, 109, 17708–17714. (41) Sockalingam, S.; Chowdhury, S. C.; Gillespie, J. W.; Keefe, M. Recent Advances in Modeling and Experiments of Kevlar Ballistic Fibrils, Fibers, Yarns and Flexible Woven Textile Fabrics – a Review. Text. Res. J. 2017, 87, 984–1010. (42) Grujicic, M.; Glomski, P. S.; He, T.; Arakere, G.; Bell, W. C.; Cheeseman, B. A. Material Modeling and Ballistic-Resistance Analysis of Armor-Grade Composites Reinforced with High-Performance Fibers. J. Mater. Eng. Perform. 2009, 18, 1169–1182. (43) Depner, M.; Schürmann, B. L. Computer Simulation of Aromatic Polyesters Including Molecular Dynamics. Polymer (Guildf). 1992, 33, 398–404. 109 (44) Brauckmann, J. O.; Zolfaghari, P.; Verhoef, R.; Klop, E. A.; de Wijs, G. A.; Kentgens, A. P. M. M. Structural Studies of Polyaramid Fibers: Solid-State NMR and First-Principles Modeling. Macromolecules 2016, 49, 5548–5560. (45) Avanzini, L.; Brambilla, L.; Marano, C.; Milani, A. Strain-Dependent Vibrational Spectra and Elastic Modulus of Poly(p-Phenylene Terephtalamide) from First-Principles Calculations. Polymer (Guildf). 2017, 116, 133–142. (46) Shimamura, K.; Misawa, M.; Ohmura, S.; Shimojo, F.; Kalia, R. K.; Nakano, A.; Vashishta, P. Crystalline Anisotropy of Shock-Induced Phenomena: Omni-Directional Multiscale Shock Technique. Appl. Phys. Lett. 2016, 108, 071901. (47) Sorin, B.; Fried, L. E. Shock Waves Science and Technology Library, Vol. 6; Zhang, F., Ed.; Springer Berlin Heidelberg: Berlin, Heidelberg, 2012. (48) Manaa, M. R.; Gee, R. H.; Fried, L. E. Internal Rotation of Amino and Nitro Groups in TATB: MP2 Versus DFT (B3LYP). J. Phys. Chem. A 2002, 106, 8806–8810. (49) Manaa, M. R. Shear-Induced Metallization of Triamino-Trinitrobenzene Crystals. Appl. Phys. Lett. 2003, 83, 1352–1354. (50) Dlott, D. D.; Fayer, M. D. Shocked Molecular Solids: Vibrational up Pumping, Defect Hot Spot Formation, and the Onset of Chemistry. J. Chem. Phys. 1990, 92, 3798–3812. (51) White, C. T.; Sinnott, S. B.; Mintmire, J. W.; Brenner, D. W.; Robertson, D. H. Chemistry and Phase Transitions from Hypervelocity Impacts. Int. J. Quantum Chem. 1994, 52, 129– 137. (52) Wu, C. J.; Fried, L. E.; Yang, L. H.; Goldman, N.; Bastea, S. Catalytic Behaviour of Dense Hot Water. Nat. Chem. 2009, 1, 57–62. (53) Zeman, S. Thermal Stabilities of Polynitroaromatic Compounds and Their Derivatives. 110 Thermochim. Acta 1979, 31, 269–283. (54) Dobratz, B. M. LLNL Explosives Handbook: Properties of Chemical Explosives and Explosives and Explosive Simulants; United States, 1981. (55) Kennedy, J. E. Second-Harmonic Generation and the Shock Sensitivity of TATB. In AIP Conference Proceedings; AIP, 2000; Vol. 505, pp 711–714. (56) Shimamura, K.; Misawa, M.; Li, Y.; Kalia, R. K.; Nakano, A.; Shimojo, F.; Vashishta, P. A Crossover in Anisotropic Nanomechanochemistry of van Der Waals Crystals. Appl. Phys. Lett. 2015, 107, 231903. (57) Tokura, Y.; Kawasaki, M.; Nagaosa, N. Emergent Functions of Quantum Materials. Nat. Phys. 2017, 13, 1056–1068. (58) Muralt, P.; Polcawich, R. G.; Trolier-McKinstry, S. Piezoelectric Thin Films for Sensors, Actuators, and Energy Harvesting. MRS Bull. 2009, 34, 658–664. (59) Muralt, P.; Polcawich, R. G.; Trolier-McKinstry, S. Piezoelectric Thin Films for Sensors, Actuators, and Energy Harvesting. MRS Bull. 2009, 34, 658–664. (60) Setter, N.; Damjanovic, D.; Eng, L.; Fox, G.; Gevorgian, S.; Hong, S.; Kingon, A.; Kohlstedt, H.; Park, N. Y.; Stephenson, G. B.; et al. Ferroelectric Thin Films: Review of Materials, Properties, and Applications. J. Appl. Phys. 2006, 100, 051606. (61) Müller, J.; Böscke, T. S.; Schröder, U.; Mueller, S.; Bräuhaus, D.; Böttger, U.; Frey, L.; Mikolajick, T. Ferroelectricity in Simple Binary ZrO 2 and HfO 2. Nano Lett. 2012, 12, 4318–4323. (62) Fichtner, S.; Wolff, N.; Lofink, F.; Kienle, L.; Wagner, B. AlScN: A III-V Semiconductor Based Ferroelectric. J. Appl. Phys. 2019, 125, 114103. (63) Müller, J.; Polakowski, P.; Mueller, S.; Mikolajick, T. Ferroelectric Hafnium Oxide Based 111 Materials and Devices: Assessment of Current Status and Future Prospects. ECS J. Solid State Sci. Technol. 2015, 4, N30–N35. (64) Umari, P.; Pasquarello, A. Ab Initio Molecular Dynamics in a Finite Homogeneous Electric Field. Phys. Rev. Lett. 2002, 89, 157602. (65) Meyer, B.; Vanderbilt, D. Ab Initio Study of BaTiO 3 and PbTiO 3 Surfaces in External Electric Fields. Phys. Rev. B 2001, 63, 205426. (66) King-Smith, R. D.; Vanderbilt, D.; King-Smith, R. D.; Vanderbilt, D. Theory of Polarization of Crystalline Solids. Phys. Rev. B 1993, 47, 1651–1654. (67) Hornik, K.; Stinchcombe, M.; White, H. Multilayer Feedforward Networks Are Universal Approximators. Neural Networks 1989, 2, 359–366. (68) Raissi, M.; Karniadakis, G. E. Hidden Physics Models: Machine Learning of Nonlinear Partial Differential Equations. J. Comput. Phys. 2018, 357, 125–141. (69) Raissi, M.; Yazdani, A.; Karniadakis, G. E. Hidden Fluid Mechanics: Learning Velocity and Pressure Fields from Flow Visualizations. J. Comput. Phys. 2018, 367, 1026–1030. (70) Owhadi, H. Bayesian Numerical Homogenization. Multiscale Model. Simul. 2015, 13, 812–828. (71) Güneş Baydin, A.; Pearlmutter, B. A.; Andreyevich Radul, A.; Mark Siskind, J. Automatic Differentiation in Machine Learning: A Survey. J. Mach. Learn. Res. 2018, 18, 1–43. (72) Rahman, A. Correlations in the Motion of Atoms in Liquid Argon. Phys. Rev. 1964, 136, A405–A411. (73) Alder, B. J.; Wainwright, T. E. Studies in Molecular Dynamics. I. General Method. J. Chem. Phys. 1959, 31, 459–466. 112 (74) Gibson, J. B.; Goland, A. N.; Milgram, M.; Vineyard, G. H. Dynamics of Radiation Damage. Phys. Rev. 1960, 120, 1229–1253. (75) Van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396–9409. (76) Liu, L.; Liu, Y.; Zybin, S. V.; Sun, H.; Goddard, W. A. ReaxFF- l g: Correction of the ReaxFF Reactive Force Field for London Dispersion, with Applications to the Equations of State for Energetic Materials. J. Phys. Chem. A 2011, 115, 11016–11022. (77) Nomura, K.; Kalia, R. K.; Nakano, A.; Vashishta, P. A Scalable Parallel Algorithm for Large-Scale Reactive Force-Field Molecular Dynamics Simulations. Comput. Phys. Commun. 2008, 178, 73–87. (78) Bloch, F. �ber Die Quantenmechanik Der Elektronen in Kristallgittern. Zeitschrift f �r Phys. 1929, 52, 555–600. (79) Phillips, J. C.; Kleinman, L. New Method for Calculating Wave Functions in Crystals and Molecules. Phys. Rev. 1959, 116, 287–294. (80) Kresse, G.; Hafner, J. Norm-Conserving and Ultrasoft Pseudopotentials for First-Row and Transition Elements. J. Phys. Condens. Matter 1994, 6, 8245–8257. (81) Troullier, N.; Martins, J. L. Efficient Pseudopotentials for Plane-Wave Calculations. II. Operators for Fast Iterative Diagonalization. Phys. Rev. B 1991, 43, 8861–8869. (82) Vanderbilt, D. Soft Self-Consistent Pseudopotentials in a Generalized Eigenvalue Formalism. Phys. Rev. B 1990, 41, 7892–7895. (83) Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953–17979. (84) Kleinman, L.; Bylander, D. M. Efficacious Form for Model Pseudopotentials. Phys. Rev. Lett. 1982, 48, 1425–1428. 113 (85) Perdew, J. P.; Zunger, A. Self-Interaction Correction to Density-Functional Approximations for Many-Electron Systems. Phys. Rev. B 1981, 23, 5048–5079. (86) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (87) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133–A1138. (88) Car, R.; Parrinello, M. Unified Approach for Molecular Dynamics and Density-Functional Theory. Phys. Rev. Lett. 1985, 55, 2471–2474. (89) Shimojo, F.; Fukushima, S.; Kumazoe, H.; Misawa, M.; Ohmura, S.; Rajak, P.; Shimamura, K.; Bassman, L.; Tiwari, S.; Kalia, R. K.; et al. QXMD: An Open-Source Program for Nonadiabatic Quantum Molecular Dynamics. SoftwareX 2019, 10, 100307. (90) Tung, I.-C.; Krishnamoorthy, A.; Sadasivam, S.; Zhou, H.; Zhang, Q.; Seyler, K. L.; Clark, G.; Mannebach, E. M.; Nyby, C.; Ernst, F.; et al. Anisotropic Structural Dynamics of Monolayer Crystals Revealed by Femtosecond Surface X-Ray Scattering. Nat. Photonics 2019, 13, 425–430. (91) Doltsinis, N. L.; Marx, D. Nonadiabatic Car-Parrinello Molecular Dynamics. Phys. Rev. Lett. 2002, 88, 166402. (92) Duncan, W. R.; Craig, C. F.; Prezhdo, O. V. Time-Domain Ab Initio Study of Charge Relaxation and Recombination in Dye-Sensitized TiO2. J. Am. Chem. Soc. 2007, 129, 8528–8543. (93) Ohmura, S.; Koga, S.; Akai, I.; Shimojo, F.; Kalia, R. K.; Nakano, A.; Vashishta, P. Atomistic Mechanisms of Rapid Energy Transport in Light-Harvesting Molecules. Appl. Phys. Lett. 2011, 98, 113302. 114 (94) Mou, W.; Ohmura, S.; Hattori, S.; Nomura, K.; Shimojo, F.; Nakano, A. Enhanced Charge Transfer by Phenyl Groups at a Rubrene/C60 Interface. J Chem Phys 2012, 136, 184705. (95) Tully, J. C. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990, 93, 1061–1071. (96) Schmidt, J. R.; Parandekar, P. V.; Tully, J. C. Mixed Quantum-Classical Equilibrium: Surface Hopping. J. Chem. Phys. 2008, 129, 044104. (97) Gross, E. K. U.; Kohn, W. Local Density-Functional Theory of Frequency-Dependent Linear Response. Phys. Rev. Lett. 1985, 55, 2850–2852. (98) Vasiliev, I.; Öğüt, S.; Chelikowsky, J. R. Ab Initio Excitation Spectra and Collective Electronic Response in Atoms and Clusters. Phys. Rev. Lett. 1999, 82, 1919–1922. (99) Yam, C.; Yokojima, S.; Chen, G. Linear-Scaling Time-Dependent Density-Functional Theory. Phys. Rev. B 2003, 68, 153105. (100) Meng, S.; Kaxiras, E. Real-Time, Local Basis-Set Implementation of Time-Dependent Density Functional Theory for Excited State Dynamics Simulations. J. Chem. Phys. 2008, 129, 054110. (101) Soto-Verdugo, V.; Metiu, H.; Gwinn, E. The Properties of Small Ag Clusters Bound to DNA Bases. J. Chem. Phys. 2010, 132, 195102. (102) Tapavicza, E.; Tavernelli, I.; Rothlisberger, U. Trajectory Surface Hopping within Linear Response Time-Dependent Density-Functional Theory. Phys. Rev. Lett. 2007, 98, 023001. (103) Walter, M.; Häkkinen, H.; Lehtovaara, L.; Puska, M.; Enkovaara, J.; Rostgaard, C.; Mortensen, J. J. Time-Dependent Density-Functional Theory in the Projector Augmented- Wave Method. J. Chem. Phys. 2008, 128, 244101. (104) Casida, M. E.; Huix-Rotllant, M. Progress in Time-Dependent Density-Functional Theory. 115 Annu. Rev. Phys. Chem. 2012, 63, 287–323. (105) Casida, M. E. Recent Advances in Density Functional Methods; Chong, D. P., Ed.; Recent Advances in Computational Chemistry; World Scientific Publishing Co. Pte. Ltd., 1995. (106) Rohlfing, M.; Louie, S. G. Electron-Hole Excitations in Semiconductors and Insulators. Phys. Rev. Lett. 1998, 81, 2312–2315. (107) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid Functionals Based on a Screened Coulomb Potential. J. Chem. Phys. 2003, 118, 8207–8215. (108) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. A Long-Range-Corrected Time-Dependent Density Functional Theory. J. Chem. Phys. 2004, 120, 8425–8433. (109) Hybertsen, M. S.; Louie, S. G. First-Principles Theory of Quasiparticles: Calculation of Band Gaps in Semiconductors and Insulators. Phys. Rev. Lett. 1985, 55, 1418–1421. (110) Dreuw, A.; Head-Gordon, M. Failure of Time-Dependent Density Functional Theory for Long-Range Charge-Transfer Excited States: The Zincbacteriochlorin−Bacteriochlorin and Bacteriochlorophyll−Spheroidene Complexes. J. Am. Chem. Soc. 2004, 126, 4007– 4016. (111) Zhang, X.; Li, Z.; Lu, G. A Non-Self-Consistent Range-Separated Time-Dependent Density Functional Approach for Large-Scale Simulations. J. Phys. Condens. Matter 2012, 24, 205801. (112) Martyna, G. J.; Tuckerman, M. E. A Reciprocal Space Based Method for Treating Long Range Interactions in Ab Initio and Force-Field-Based Calculations in Clusters. J. Chem. Phys. 1999, 110, 2810–2821. (113) Mou, W.; Hattori, S.; Rajak, P.; Shimojo, F.; Nakano, A. Nanoscopic Mechanisms of Singlet Fission in Amorphous Molecular Solid. Appl. Phys. Lett. 2013, 102, 173301. 116 (114) Reed, E. J.; Fried, L. E.; Joannopoulos, J. D. A Method for Tractable Dynamical Studies of Single and Double Shock Compression. Phys. Rev. Lett. 2003, 90, 235503. (115) Resta, R.; Vanderbilt, D. Theory of Polarization: A Modern Approach. In Physics of Ferroelectrics: A Modern Perspective; Springer Berlin Heidelberg: Berlin, Heidelberg, 2007; pp 31–68. (116) Raffaele Resta. Manifestations of Berry’s Phase in Molecules and Condensed Matter. J. Phys. Condens. Matter 2000, 12, R107. (117) Schumacher, M.; Weber, H.; Jóvári, P.; Tsuchiya, Y.; Youngs, T. G. A.; Kaban, I.; Mazzarello, R. Structural, Electronic and Kinetic Properties of the Phase-Change Material Ge2Sb2Te5 in the Liquid State. Sci. Rep. 2016, 6, 27434. (118) Waldecker, L.; Miller, T. A.; Rudé, M.; Bertoni, R.; Osmond, J.; Pruneri, V.; Simpson, R. E.; Ernstorfer, R.; Wall, S. Time-Domain Separation of Optical Properties from Structural Transitions in Resonantly Bonded Materials. Nat. Mater. 2015, 14, 991–995. (119) Kolobov, A. V.; Krbal, M.; Fons, P.; Tominaga, J.; Uruga, T. Distortion-Triggered Loss of Long-Range Order in Solids with Bonding Energy Hierarchy. Nat. Chem. 2011, 3, 311– 316. (120) Rousse, A.; Rischel, C.; Fourmaux, S.; Uschmann, I.; Sebban, S.; Grillon, G.; Balcou, P.; Förster, E.; Geindre, J. P.; Audebert, P.; et al. Non-Thermal Melting in Semiconductors Measured at Femtosecond Resolution. Nature 2001, 410, 65–68. (121) Kolobov, A. V; Tominaga, J.; Fons, P. Springer Handbook of Electronic and Photonic Materials. Springer Handb. Electron. Photonic Mater. 2017, 1149–1161. (122) Liu, B.; Song, Z.; Feng, S.; Chen, B. Characteristics of Chalcogenide Nonvolatile Memory Nano-Cell-Element Based on Sb2Te3 Material. Microelectron. Eng. 2005, 82, 168–174. 117 (123) Zijlstra, E. S.; Kalitsov, A.; Zier, T.; Garcia, M. E. Fractional Diffusion in Silicon. Adv. Mater. 2013, 25, 5605–5608. (124) Jeyasingh, R.; Fong, S. W.; Lee, J.; Li, Z.; Chang, K.-W.; Mantegazza, D.; Asheghi, M.; Goodson, K. E.; Wong, H.-S. P. Ultrafast Characterization of Phase-Change Material Crystallization Properties in the Melt-Quenched Amorphous Phase. Nano Lett. 2014, 14, 3419–3426. (125) Poborchii, V. V.; Kolobov, A. V.; Tanaka, K. Photomelting of Selenium at Low Temperature. Appl. Phys. Lett. 1999, 74, 215–217. (126) Fons, P.; Osawa, H.; Kolobov, A. V.; Fukaya, T.; Suzuki, M.; Uruga, T.; Kawamura, N.; Tanida, H.; Tominaga, J. Photoassisted Amorphization of the Phase-Change Memory Alloy Ge2Sb2Te 5. Phys. Rev. B 2010, 82, 041203. (127) Saeta, P.; Wang, J.-K.; Siegal, Y.; Bloembergen, N.; Mazur, E. Ultrafast Electronic Disordering during Femtosecond Laser Melting of GaAs. Phys. Rev. Lett. 1991, 67, 1023– 1026. (128) Stampfli, P.; Bennemann, K. H. Dynamical Theory of the Laser-Induced Lattice Instability of Silicon. Phys. Rev. B 1992, 46, 10686–10692. (129) CRC Handbook of Chemistry and Physics, 95th ed.; Haynes, W. M., Ed.; CRC Press: Boca Raton, 2014. (130) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511–519. (131) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Explicit Reversible Integrators for Extended Systems Dynamics. Mol. Phys. 1996, 87, 1117–1157. (132) Shim, V. P. W.; Lim, C. T.; Foo, K. J. Dynamic Mechanical Properties of Fabric Armour. 118 Int. J. Impact Eng. 2001, 25, 1–15. (133) Penn, L.; Larsen, F. Physicochemical Properties of Kevlar 49 Fiber. J. Appl. Polym. Sci. 1979, 23, 59–73. (134) Deteresa, S. J.; Allen, S. R.; Farris, R. J.; Porter, R. S. Compressive and Torsional Behaviour of Kevlar 49 Fibre. J. Mater. Sci. 1984, 19, 57–72. (135) Qi, L.; Sinnott, S. B. Polymerization via Cluster−Solid Surface Impacts: Molecular Dynamics Simulations. J. Phys. Chem. B 1997, 101, 6883–6890. (136) Patterson, J. E.; Lagutchev, A.; Huang, W.; Dlott, D. D. Ultrafast Dynamics of Shock Compression of Molecular Monolayers. Phys. Rev. Lett. 2005, 94, 015501. (137) Kadau, K. Microscopic View of Structural Phase Transitions Induced by Shock Waves. Science (80-. ). 2002, 296, 1681–1684. (138) Kushner, H. J. A New Method of Locating the Maximum Point of an Arbitrary Multipeak Curve in the Presence of Noise. J. Basic Eng. 1964, 86, 97–106. (139) Rasmussen, C. E.; Williams, C. K. I. Gaussian Processes for Machine Learning.; MIT Press Cambridge, MA, 2004; Vol. 14. (140) Sun, J.; Ruzsinszky, A.; Perdew, J. Strongly Constrained and Appropriately Normed Semilocal Density Functional. Phys. Rev. Lett. 2015, 115, 036402. (141) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long- Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787–1799. (142) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (143) Sabatini, R.; Gorni, T.; de Gironcoli, S. Nonlocal van Der Waals Density Functional Made 119 Simple and Efficient. Phys. Rev. B 2013, 87, 041108. (144) Keten, S.; Xu, Z.; Ihle, B.; Buehler, M. J. Nanoconfinement Controls Stiffness, Strength and Mechanical Toughness of β-Sheet Crystals in Silk. Nat. Mater. 2010, 9, 359–367. (145) Strachan, A.; van Duin, A. C. T.; Chakraborty, D.; Dasgupta, S.; Goddard, W. A. Shock Waves in High-Energy Materials: The Initial Chemical Events in Nitramine RDX. Phys. Rev. Lett. 2003, 91, 098301. (146) Shao, Y.; Gan, Z.; Epifanovsky, E.; Gilbert, A. T. B.; Wormit, M.; Kussmann, J.; Lange, A. W.; Behn, A.; Deng, J.; Feng, X.; et al. Advances in Molecular Quantum Chemistry Contained in the Q-Chem 4 Program Package. Mol. Phys. 2015, 113, 184–215. (147) Yeh, W. Y.; Young, R. J. Molecular Deformation Processes in Aromatic High Modulus Polymer Fibres. Polymer (Guildf). 1999, 40, 857–870. (148) Zhang, Y.; Guo, G.; Refson, K.; Zhao, Y. Finite-Size Effect at Both High and Low Temperatures in Molecular Dynamics Calculations of the Self-Diffusion Coefficient and Viscosity of Liquid Silica. J. Phys. Condens. Matter 2004, 16, 9127–9135. (149) Heyes, D. M.; Cass, M. J.; Powles, J. G.; Evans, W. A. B. Self-Diffusion Coefficient of the Hard-Sphere Fluid: System Size Dependence and Empirical Correlations. J. Phys. Chem. B 2007, 111, 1455–1464. (150) Schnell, S. K.; Liu, X.; Simon, J.-M.; Bardow, A.; Bedeaux, D.; Vlugt, T. J. H.; Kjelstrup, S. Calculating Thermodynamic Properties from Fluctuations at Small Scales. J. Phys. Chem. B 2011, 115, 10911–10918. (151) Krüger, P.; Schnell, S. K.; Bedeaux, D.; Kjelstrup, S.; Vlugt, T. J. H.; Simon, J.-M. Kirkwood–Buff Integrals for Finite Volumes. J. Phys. Chem. Lett. 2013, 4, 235–238. (152) Shen, Y.; Jester, S. B.; Qi, T.; Reed, E. J. Nanosecond Homogeneous Nucleation and 120 Crystal Growth in Shock-Compressed SiO2. Nat. Mater. 2016, 15, 60–65. (153) Manaa, M. R.; Reed, E. J.; Fried, L. E.; Goldman, N. Nitrogen-Rich Heterocycles as Reactivity Retardants in Shocked Insensitive Explosives. J. Am. Chem. Soc. 2009, 131, 5483–5487. (154) Östmark, H. Shock Induced Sub-Detonation Chemical Reactions in 1,3,5-Triamino-2,4,6- Trinitrobenzene. AIP Conf. Proc. 1996, 370, 871–874. (155) Wu, C. J.; Yang, L. H.; Fried, L. E.; Quenneville, J.; Martinez, T. J. Electronic Structure of Solid 1,3,5-Triamino-2,4,6-Trinitrobenzene under Uniaxial Compression: Possible Role of Pressure-Induced Metallization in Energetic Materials. Phys. Rev. B 2003, 67, 235101. (156) Farber, M.; Srivastava, R. D. Thermal Decomposition of 1,3,5-Triamino-2,4,6- Trinitrobenzene. Combust. Flame 1981, 42, 165–171. (157) Sharma, J.; Garrett, W. L.; Owens, F. J.; Vogel, V. L. X-Ray Photoelectron Study of Electronic Structure, Ultraviolet, and Isothermal Decomposition of 1,3,5-Triamino-2,4,6- Trinitrobenzene. J. Phys. Chem. 1982, 86, 1657–1661. (158) Cady, H. H.; Larson, A. C. The Crystal Structure of 1,3,5-Triamino-2,4,6-Trinitrobenzene. Acta Crystallogr. 1965, 18, 485–496. (159) Sharma, J.; Forbes, J. W.; Coffey, C. S.; Liddiard, T. P. The Physical and Chemical Nature of Sensitization Centers Left from Hot Spots Caused in Triaminotrinitrobenzene by Shock or Impact. J. Phys. Chem. 1987, 91, 5139–5144. (160) Makashir, P. S.; Kurian, E. M. Spectroscopic and Thermal Studies on the Decomposition of l,3,5-Triamino-2,4,6-Trinitrobenzene (TATB). J. Therm. Anal. 1996, 46, 225–236. (161) Zhang, L.; Zybin, S. V; van Duin, A. C. T.; Dasgupta, S.; Goddard, W. A.; Kober, E. M. Carbon Cluster Formation during Thermal Decomposition of Octahydro-1,3,5,7- 121 Tetranitro-1,3,5,7-Tetrazocine and 1,3,5-Triamino-2,4,6-Trinitrobenzene High Explosives from ReaxFF Reactive Molecular Dynamics Simulations. J. Phys. Chem. A 2009, 113, 10619–10640. (162) Quenneville, J.; Germann, T. C.; Thompson, A. P.; Kober, E. M. Molecular Dynamics Studies of Thermal Induced Chemistry in TATB. AIP Conf. Proc. 2007, 955, 451–454. (163) Li, Y.; Kalia, R. K.; Nakano, A.; Nomura, K. I.; Vashishta, P. Multistage Reaction Pathways in Detonating High Explosives. Appl. Phys. Lett. 2014, 105. (164) Chang, J.; Lian, P.; Wei, D.-Q.; Chen, X.-R.; Zhang, Q.-M.; Gong, Z.-Z. Thermal Decomposition of the Solid Phase of Nitromethane: Ab Initio Molecular Dynamics Simulations. Phys. Rev. Lett. 2010, 105, 188302. (165) Schweigert, I. V. Ab Initio Molecular Dynamics of High-Temperature Unimolecular Dissociation of Gas-Phase RDX and Its Dissociation Products. J. Phys. Chem. A 2015, 119, 2747–2759. (166) Ge, N.-N.; Wei, Y.-K.; Ji, G.-F.; Chen, X.-R.; Zhao, F.; Wei, D.-Q. Initial Decomposition of the Condensed-Phase β-HMX under Shock Waves: Molecular Dynamics Simulations. J. Phys. Chem. B 2012, 116, 13696–13704. (167) Britt, A. D.; Moniz, W. B.; Chingas, G. C.; Moore, D. W.; Heller, C. A.; Ko, C. L. Free Radicales of TATB. Propellants, Explos. Pyrotech. 1981, 6, 94–95. (168) Ornellas, D. L. Calorimetric Determinations of the Heat and Products of Detonation for Explosives: October 1961 to April 1982; 1982. (169) He, Z.-H.; Chen, J.; Wu, Q. Initial Decomposition of Condensed-Phase 1,3,5-Triamino- 2,4,6-Trinitrobenzene under Shock Loading. J. Phys. Chem. C 2017, 121, 8227–8235. (170) Stavrou, E.; Crowhurst, J. C.; Zaug, J. M.; Armstrong, M. R.; Fried, L. E.; Radousky, B. 122 No Title. Present. APS March Meet. 2016, Balt. 2016. (171) Akiyama, M.; Kamohara, T.; Kano, K.; Teshigahara, A.; Takeuchi, Y.; Kawahara, N. Enhancement of Piezoelectric Response in Scandium Aluminum Nitride Alloy Thin Films Prepared by Dual Reactive Cosputtering. Adv. Mater. 2009, 21, 593–596. (172) Farrer, N.; Bellaiche, L. Properties of Hexagonal ScN versus Wurtzite GaN and InN. Phys. Rev. B - Condens. Matter Mater. Phys. 2002, 66, 2012031–2012034. (173) Tasnádi, F.; Alling, B.; Höglund, C.; Wingqvist, G.; Birch, J.; Hultman, L.; Abrikosov, I. A. Origin of the Anomalous Piezoelectric Response in Wurtzite ScxAl1-XN Alloys. Phys. Rev. Lett. 2010, 104, 1–4. (174) Feynman, R. P. The Principle of Least Action in Quantum Mechanics., Princeton University, 1942. (175) Feynman 1918-1988, R. P. (Richard P. The Feynman Lectures on Physics; Reading, Mass. : Addison-Wesley Pub. Co., c1963-1965. (176) Raissi, M.; Perdikaris, P.; Karniadakis, G. E. Physics-Informed Neural Networks: A Deep Learning Framework for Solving Forward and Inverse Problems Involving Nonlinear Partial Differential Equations. J. Comput. Phys. 2019, 378, 686–707.
Abstract (if available)
Abstract
Quantum materials have emerged as a material of choice for new computing paradigm and information storage. Emergent properties of quantum materials are generally studied under non-equilibrium conditions such as photoexcitation and high electric field. Modeling far from equilibrium condition for quantum materials is a new frontier in molecular modeling. Quantum material’s exotic properties arise from intrinsically quantum phenomena, and thus requiring a quantum treatment of non-equilibrium processes during modeling. This thesis explores the application of quantum molecular dynamics (QMD) under three different non-equilibrium conditions for various materials. ❧ First, we perform nonadiabatic quantum molecular dynamics (NAQMD) simulations to study photoexcitation-driven non-thermal amorphization of phase-change materials, GeTe and Sb₂Te₃. Phase change materials are of high interest for low-power, high-throughput storage devices for next generation neuromorphic computing technologies. Their usage is based on the contrasting properties of their amorphous and crystalline phases, which can be switched in the nanosecond timescale. Here, we show that the energy required for, and the time associated with, the amorphization process can be further reduced by using a non-thermal photoexcitation path. In our NAQMD simulations for GeTe and Sb₂Te₃, non-thermal pathway was characterized by an instantaneous charge transfer from Te-p orbitals to Ge-p/Sb-p orbitals upon excitation. Time evolution of excited state shows an increasing bonding characteristic for M-M (M = Ge, Sb) and decreasing M-Te interaction, leading to the formation of wrong bonds. These results highlight an energy-efficient and ultrafast amorphization pathway that could be used to enhance the performance of phase change materials based opto-electronic devices. ❧ Second, we model shock in aramid fibers with quantum and classical molecular dynamics (MD) simulations using multiscale shock theory (MSST). Aramid fibers are widely used as lightweight, shock-resistance, reinforcing materials. However, their intrinsic shock behaviors are not well understood. We perform ab-initio QMD simulations of shock on poly(p-phenylene terephthalamide) (PPTA), which reveals stress release mechanisms based on structural phase transformation (SPT) and generation of planar amorphous region. SPT is triggered by [100] shock-induced coplanarity of phenylene groups and rearrangement of sheet stacking, leading to a previously unknown monoclinic phase. Planar amorphous region is generated by [010] shock-induced scission of hydrogen bonds leading to disruption of polymer sheets, and trans-to-cis conformational change of polymer chains. In contrast to the latter, the former mechanism preserves hydrogen bonding and cohesiveness of polymer chains in the identified novel crystalline phase preserving the strength of PPTA. The interplay between hydrogen bond preserving (SPT) and non-preserving (planar amorphous) shock release mechanisms are critical to understanding the shock performance of aramid fibers. To understand the size effect, we also perform million-atom MD simulation of Kevlar using reactive forcefield. In [100] direction, we observed the structural phase transformation similar to that in QMD. On the contrary, in [010] direction, we found the formation of paracrystalline Kevlar for large systems, which was impossible to study by QMD. ❧ Third, we model non-equilibrium shock simulation to understand the detonation process in 2,4,6-triamino-1,3,5-trinitrobenzene (TATB). Detonation processes probed with atomistic details have remained elusive due to highly complex reactions in heterogeneous shock structures. Here, we provide atomistic details of the initial reaction pathways during shock-induced decomposition of TATB crystal using large reactive molecular dynamics simulations based on reactive force fields. Simulation results reveal the existence of three competing intermolecular pathways for the formation of N₂. We also observe the formation of large nitrogen- and oxygen-rich carbon aggregates, which delays the release of final reaction products. ❧ Finally, we use adiabatic QMD to perform computational synthesis of Nitride ferroelectrics, Al₁₋ₓScₓN (x = 0.27 – 0.43). We use a Berry-phase approach in conjunction with QMD to study electric field-induced switching of polarization in Nitride ferroelectrics. The computed electric field value for polarization switching in Al₀.₇₅Sc₀.₂₅N (600 MV/m) is in close agreement with the experimental value (400 MV/m). Furthermore, our simulations correctly predict the lowering of switching field with increasing Sc concentration. QMD also determines the switching pathway and effect of Sc doping in Al₁₋ₓScₓN. ❧ Apart from all the traditional MD techniques, machine learning (ML) has recently grown in popularity. In particular, ML has been used to study the ground state of materials. Here, we use ML to extract both ground and excited states of the hydrogen atom from the Schrödinger equation. The ML technique used to learn Schrödinger equation provides energies and wavefunctions of ground and excited states, which are in excellent agreement with the analytical results. Remarkably, the combined ML and self-consistent methods provide exact wave functions without Gram-Schmidt orthogonalization. ❧ With ever-increasing knowledge of materials, we need better tools to perform high-throughput screening of material and model fundamental and exotic properties occurring in the materials with high fidelity. In the present era of quantum materials where non-equilibrium processes are prevalent, this thesis demonstrates a roadmap to study non-equilibrium process using quantum molecular dynamics and machine learning.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Reactive and quantum molecular dynamics study of materials: from oxidation and amorphization to interfacial charge transfer
PDF
Controlling electronic properties of two-dimensional quantum materials: simulation at the nexus of the classical and quantum computing eras
PDF
Multi-scale quantum dynamics and machine learning for next generation energy applications
PDF
Large-scale molecular dynamics simulations of nano-structured materials
PDF
High throughput computational framework for synthesis and accelerated discovery of dielectric polymer materials using polarizable reactive molecular dynamics and graph neural networks
PDF
Atomistic modeling of the mechanical properties of metallic glasses, nanoglasses, and their composites
PDF
Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations
PDF
Photoexcitation and nonradiative relaxation in molecular systems: methodology, optical properties, and dynamics
PDF
Simulation and machine learning at exascale
PDF
Molecular dynamics simulations of nanostructures
PDF
Modeling of bicontinuous nanoporous materials and investigation of their mechanical and hydraulic properties
PDF
Out-of-equilibrium dynamics of inhomogeneous quantum systems
PDF
Shock sensitivity of energetic material and nanometric damage mechanisms in silica glass
PDF
Dynamics of wing cracks and nanoscale damage in silica glass
PDF
Adiabatic and non-adiabatic molecular dynamics in nanoscale systems: theory and applications
PDF
Shock-induced poration, cholesterol flip-flop and small interfering RNA transfection in a phospholipid membrane: multimillion atom, microsecond molecular dynamics simulations
PDF
Theoretical modeling of nanoscale systems: applications and method development
PDF
Molecular dynamics studies of protein aggregation in unbounded and confined media
PDF
Physics informed neural networks and electrostrictive cavitation in water
PDF
Phase change heterostructures for electronic and photonic applications
Asset Metadata
Creator
Tiwari, Subodh Chandra
(author)
Core Title
Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Materials Science
Publication Date
05/12/2020
Defense Date
03/19/2020
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
excited state dynamics,ferroelectric switching,Kevlar,machine learning,non-equilibrium process,OAI-PMH Harvest,phase change materials,quantum molecular dynamics,reactive molecular dynamics,shock simulation
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Vashishta, Priya (
committee chair
), Branicio, Paulo (
committee member
), Nakano, Aiichiro (
committee member
)
Creator Email
sctiwari@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c89-303175
Unique identifier
UC11665835
Identifier
etd-TiwariSubo-8487.pdf (filename),usctheses-c89-303175 (legacy record id)
Legacy Identifier
etd-TiwariSubo-8487.pdf
Dmrecord
303175
Document Type
Dissertation
Rights
Tiwari, Subodh Chandra
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
excited state dynamics
ferroelectric switching
Kevlar
machine learning
non-equilibrium process
phase change materials
quantum molecular dynamics
reactive molecular dynamics
shock simulation