University of Southern California Dissertations and Theses

Theoretical modeling of nanoscale systems: applications and method development

(USC Thesis Other)

Theoretical modeling of nanoscale systems: applications and method development

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Theoretical Modeling of Nanoscale Systems: Applications and Method Development By Guoqing Zhou 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 (PHYSICS) May 2021 Copyright 2021 Guoqing Zhou ii Acknowledgements I would like to express my appreciation to the many people for their support throughout my education. First of all, I would like to thank my doctoral adviser Prof. Oleg Prezhdo for his continuous support during my Ph.D. journey both pedagogically and financially. Oleg is a great research adviser who taught me how to find valuable scientific problems and approach them in systematic ways. He is also an effective mentor who taught me the value of challenging oneself, overcoming weakness for perfectness. Thank you for your mentorship which has shaped me into the persistent and critical thinker that I am today. I am grateful to the professors in the CACS group in USC. They are Prof. Aiichiro Nakano, Prof. Priya Vashishta and my previous adviser Prof. Rajiv Kalia. Thanks for their support and suggestion which guided me into the field of computational physics. Prof. Kalia’s patience and motivation are invaluable to me during the hard time when I just started my research projects. I would like to thank the following people with whom I closely collaborated with. They are Sergei Tretiak (staff scientist, LANL), Benjamin Nebgen (staff scientist, LANL), Nicholas Lubber (staff scientist, LANL), Anders Niklasson (staff scientist, LANL), Dongyan Xu (Professor, CUHK), Mingsen Deng (Professor, Guizhou Education University), Pulickel Ajayan (Professor, Rice), Weibin Chu (postdoc, USC), Yucheng Xiong (postdoc, CUHK), Walter Malone (postdoc, LANL), Pankaj Rajak (postdoc, Argonne National Laboratory). Their fruitful discussion and suggestions, constructive criticisms have made this dissertation possible. I am grateful to my undergraduate professors at University of Science and Technology of China who cultivated my interest in science. Prof. Banjiao Ye is my first physics teacher, his systematic and attractive introduction courses sparked my interest in physics and led to this journey of pursuing the Ph.D. Prof. Zejun Ding’s course in computational physics was both conceptually iii and mathematically stimulating. His classes were full of open questions which challenged us to think creatively. It was his invaluable introduction which led me to computational physics in graduate school. A special thanks to the members of my qualifying and thesis committees: Professors Oleg Prezhdo, Krzysztof Pilch, Stephan Haas, Curt Wittig, and James Boedicker. I appreciate you making the time given your busy schedules. Thank you! I would like to thank my family and friends. A great thanks to my brother Guodong Zhou. He is a postdoc in CUHK and help me connect and find collaboration with researchers there. Most importantly, he introduced me into science with a range of science fictions, fundamental science magazines, which planted the seed of science in my mind. I thank my parents and friends, for their support and accompaniment over the years and for the many years to come. Finally, I am grateful to the NSF, DOE and LANL for financially supporting the research presented in this dissertation. iv Table of Contents Acknowledgements ........................................................................................................................ ii List of Tables ................................................................................................................................. vi List of Figures ............................................................................................................................... vii Abstract ......................................................................................................................................... ix Chapter 1: Introduction ................................................................................................................... 1 1.1 Classical and Born-Oppenheimer Molecular Dynamics ................................... 3 1.2 Non-Adiabatic Molecular Dynamics with Surface Hopping ............................ 5 1.3 Machine Learning in Computational Material Science .................................... 7 1.4 Thesis Organization .......................................................................................... 9 Chapter 2: Molecular Simulation of MoS2 Exfoliation ................................................................. 11 2.1 Introduction .................................................................................................... 11 2.2 Simulation Setup ............................................................................................. 13 2.3 Results and Discussion .................................................................................... 15 2.3.1 Shock Impact on Solvent .................................................................. 15 2.3.2 Shocked Induced Bubble Collapse ................................................... 17 2.3.3 Onset of exfoliation .......................................................................... 18 2.3.4 Exfoliation with repeated impact of shock wave .............................. 22 2.3.5 Experimental Results ....................................................................... 23 2.4 Conclusions .................................................................................................... 24 2.5 Method ............................................................................................................ 25 Chapter 3: A Comparative Study of Electron–Phonon Scattering in Carbon Nanotubes and Graphene Nanoribbons ............................................................................................... 29 3.1 Introduction .................................................................................................... 29 3.2 Computational Details .................................................................................... 32 3.3 Results and Discussion .................................................................................... 34 3.4 Conclusions .................................................................................................... 43 Chapter 4: Modeling Auger Processes with Non-Adiabatic Molecular Dynamics ....................... 45 4.1 Introduction .................................................................................................... 45 4.2 Computational Details .................................................................................... 48 4.3 Case Study on CdSe Quantum Dot .................................................................. 52 4.4 Conclusions .................................................................................................... 60 Chapter 5: Unsupervised Machine Learning of Non-Adiabatic Molecular Dynamics in MAPbI3 63 5.1 Introduction .................................................................................................... 63 5.2 Computational Details .................................................................................... 66 5.2.1 Non-Adiabatic Molecular Dynamics .............................................. 66 5.2.2 Mutual Information .......................................................................... 68 v 5.3 Results and Discussion .................................................................................... 69 5.4 Conclusions .................................................................................................... 77 Chapter 6: GPU-Accelerated Semi-Empirical Born Oppenheimer Molecular Dynamics using PyTorch ....................................................................................................................... 79 6.1 Introduction .................................................................................................... 79 6.2 Methods .......................................................................................................... 83 6.2.1 Semi-Empirical Quantum Chemistry Methods ................................ 83 6.2.2 SP2 algorithm .................................................................................. 88 6.2.3 Extended Lagrangian BOMD .......................................................... 90 6.3 Results and Discussion .................................................................................... 93 6.3.1 Nanostar Dendrimer ......................................................................... 93 6.3.2 Polyethylene .................................................................................... 98 6.4 Conclusions .................................................................................................. 103 Chapter 7: Machine Learning with Domain Knowledge in Quantum Chemistry......................... 106 7.1 Introduction .................................................................................................. 106 7.2 Model ............................................................................................................ 108 7.3 Training Details ............................................................................................ 111 7.4 Results and Discussion .................................................................................. 112 7.5 Conclusions .................................................................................................. 118 Chapter 8: Closing Remarks ……………................................................................................... 120 References .................................................................................................................................. 123 Appendices ................................................................................................................................. 142 Appendix A: Supplementary Information for Chapter 2 ..................................... 142 A1. Force-Field Validation by experimental measurement ..................... 142 A2. Force-Field Validation for Shock Simulation ................................... 143 A3. Additional Results ............................................................................ 144 Appendix B: Supporting Information for Chapter 4 ........................................... 147 Appendix C: Supporting Information for Chapter 5 ............................................ 149 C1. Theoretical Methodology ................................................................. 149 C2. Features from the Molecular Dynamics trajectories ......................... 150 C3. Unsupervised Machine Learning with Mutual Information .............. 152 vi List of Tables 3.1 Intra-band relaxation times for hot electrons and holes in the (6,5) Carbon nanotubes (SWCNT) and graphene nanoribbons (GNR). P-37 3.2 Electron-hole recombination times, pure-dephasing times, and average absolute nonadiabatic coupling (NAC) in the SWCNT and GNR. P-38 4.1 Averaged Coulomb matrix elements and absolute values of NAC matrix elements. P-58 4.2 Timescales for hot electron relaxion and charge recombination. P-59 5.1 Mutual information of internal coordinates and motions with NAC. P-72 5.2 Mutual information of internal coordinates with bandgap. P-75 6.1 Energy drift and fluctuation for Born-Oppenheimer Molecular Dynamics (BOMD) and extended Lagrangian BOMD. P-98 7.1 Accuracy on predicting atomization energies and atomic forces. P-116 S2.1 Contact angles calculated by MD simulations for various mass fractions of IPA in the mixture. P-143 S5.1 Static features used in the mutual information analysis. P-150 vii List of Figures 2.1 Initial configuration of the shock-induced exfoliation simulation and formation of ice VII motif. P-14 2.2 Snapshots of the collapsing cavitation bubble. P-17 2.3 Collapse of Cavitation bubble and the generation of a nanojet and a vortex. P-19 2.4 Pressure, temperature and shear stress distributions in MoS2 after the nanojet impact. P-20 2.5 Exfoliated MoS2 with solvent between the nanosheets. P-23 3.1 Simulation cells with optimized structures of (6,5) carbon nanotube (SWCNT), and graphene nanoribbon (GNR). P-32 3.2 Geometric structure, and HOMO and LUMO density distributions of the SWCNT and GNR at 0 K and 300 K. P-35 3.3 Intra-band relaxation of hot electrons and holes in the SWCNT and GNR. P-37 3.4 Nonadiabatic coupling between Kohn-Sham orbitals in SWCNT and GNR. P-40 3.5 Electron-phonon influence spectra for electron-hole recombination in SWCNT and GNR. P-42 4.1 Structure and electronic structure of the CdSe quantum dot (QD) and a schematic view of Auger process. P-54 4.2 State populations illustrating the Auger process. P-56 4.3 NAC between valence and conduction band orbitals. P-57 4.4 Decay of the total electronic energy in CdSe QD. P-59 4.5 Energy transfer from electron to hole during Auger process and the phonon influence spectrum for hole relaxation. P-60 viii 5.1 Geometric and electronic structure of MAPbI3. P-66 5.2 Probability distributions of geometric and electronic quantities. P-70 5.3 Evolution of bond angles and NAC and their joint distribution. P-71 5.4 NAC and bandgap spectra compared with spectra of the key internal coordinates. P-76 6.1 Stacked histogram of the average time spent computing several key quantities. P-95 6.2 Energy drift of the nanostar system during simulations. P-97 6.3 Average time spending on computing key quantities. P-100 6.4 Average time spending on computing key quantities with the batch mode enabled. P-102 6.5 Performance boost with the batch mode compared with running in serial. P-103 7.1 Model structure scheme. P-108 7.2 Distribution of molecule size used in training, testing. P-112 7.3 Output SEQM parameter distributions from hybrid model. P-114 7.4 Mean absolute error vs molecule sizes for atomization energies and atomic forces P-118 S2.1 Side view of the system for the calculation of contact angle. P-143 S2.2 Hugoniot compression curves of water and water/IPA mixture. P-144 S2.3 Exfoliated MoS2 with solvent between the nanosheets. P-145 S2.4 Exfoliation of MoS2 samples. P-146 S4.1 Energy decay of hot electrons in CdSe QD. P-147 S4.2 Energy decay of hot holes in CdSe QD. P-148 S5.1 Geometric structure of MAPbI3 to show bonds and angles. P-152 ix Abstract Nanomaterials exhibit extraordinary thermal, mechanic, optical, electrical and magnetic properties and have wide applications in engineering, materials science and physical chemistry. The theoretical modeling of these nanoscale systems provides a fundamental understanding of various processes, such as physical and chemical synthesis, nuclei and electron dynamics, carrier and energy transfer. Based on the mechanism and timescale, the modeling of these processes requires the development and utilization of theories in classical and quantum dynamics, electronic structure. In this dissertation, two main methods, classical molecular dynamics and non-adiabatic molecular dynamics, will be presented in different scenarios with applied on studying the liquid phase exfoliation process in generating layered materials and the carrier dynamics in carbon nanotube and graphene nanoribbon. We detail the computational schemes along with accurate description for properties of interest. At the same time, we highlight their limitations and disadvantages and the contributions we make to improve their applicability. One contribution we make is to incorporate proper carrier interactions for modeling multi-particle processes with non- adiabatic molecular dynamics. Apart from this, due to the dense manifold of electronic degrees of freedom in nanoscale systems, the modeling requires efficient algorithms and software, and accurate analyzing tools. We incorporate techniques from machine learning into the theoretical modeling on nanoscale systems. We present an unsupervised learning technique for extracting implicit dependency of properties affecting nuclei and electronic dynamics and provides quantitative estimation and comparison of impact from lattice system on the electron recombination. In addition, we present a semi-empirical quantum chemistry model implemented with machine learning framework for improving efficiency and feasibility. The package is designed to take advantage of modern computer architecture graphic process units and is x benchmarked with a set of organic molecules. Finally, we present a hybrid model with the developed package and supervised learning technique: artificial neural network. The hybrid model is developed to address the interpretability and transferability problem in machine learning when applying into material science. Incorporating the domain knowledge of quantum chemistry helps alleviate these issues and at the same time, machine learning greatly improves the accuracy of this semi-empirical model. The simplicity, efficiency and accuracy from this hybrid semi-empirical model allows the applications on large molecule systems for modeling long time processes and requiring quantum chemistry level accuracy. This model can be further extended to incorporate with non-adiabatic molecular dynamics for studying processes involving fast and rapid electron dynamics and to push forward the edges of theoretical modeling on nanoscale systems. 1 Chapter 1: Introduction Nanomaterials exhibit unique structural and electrical properties with sizes ranging from 1 to 100 nm. Based on their sizes and the performance of modern computer architectures, the modeling of these systems can have a full atomistic description and provide insight understanding of the systems, help illustrate the materials’ interaction and explain mechanisms of dynamics. As such, computational modeling is important to both theoretical and experimental research. There are many cases that experiments are difficult to fully understand the processes due to the time- and length-scales, also it may be costly and time-consuming to perform experiments. In these cases, theoretical modeling can help relieve the experiment work. For instance, in drug design, theoretical modeling can help narrow down the search of candidate drug molecules with identifying typical properties, or in the design of photovoltaic materials, modeling can unveil the origin of performance bottleneck and explore the potential modification on the systems before the physical synthesis. Theoretical modeling plays a central role in material science and chemistry and connects with experiments for the study of various materials and the exploration of natural science. Two main groups of methods have been developed since several decades ago for the study of processes in nanoscale systems with timescales ranging from femtosecond to sub-microseconds. One is classical molecular dynamics (CMD), 1 which assumes the electrons on ground state, neglects the electron degrees of freedom, and approximates the interactions with empirical formulas. Because of these approximations, CMD is capable of modeling systems up to billions of atoms with current modern supercomputers and studying processes with timescale up to 10-100 nanoseconds. In CMD, the motions of atoms and molecules are determined by numerically solving Newton’s equations of motion and the interaction of particles is described by formulas or force fields like harmonic bonds, bond-angles, screened Coulomb and vdw non-bonds, etc. CMD has 2 been successfully applied to study the processes like chemical and physical synthesis, 2-3 atom diffusion, 4 heat transport, 5 phase transition, 6 etc. However, the exclusion of electron degrees of freedom restricts the CMD from modeling processes involving electron dynamics. Due to the simplification of formula terms describing the interaction between particles, CMD usually requires force field development to tune the properties of interest not only for studying different systems, but also for modeling different processes. One more accurate but computationally demanding method is Born-Oppenheimer Molecular Dynamics (BOMD). 7 In BOMD, it assumes that the fast- evolving electrons will quickly adapt the slow motion of nuclei and stay in ground state, as such their motions can be treated separately. With this adiabatic approximation, the motion of nuclei is governed by the ground state potential energy surface from first principle calculation. Because of the cost of first principle calculation, BOMD is limited to use on small systems with up to hundreds of atoms for processes with timescales up to 10-100 picoseconds. When modeling processes involving electron dynamics, the adiabatic approximation breaks down. Non-adiabatic molecular dynamics (NAMD) will be preferred in these cases where electronic transition may happen due to strong electron-nuclear coupling or nonadiabatic coupling (NAC). There are several methods developed, 8-9 and among them Tully’s surface hopping (SH) 10 is the most popular one for its efficiency and simplicity. Many extensions of SH have been developed to address problems like decoherance, 11 detail balance, 12 and modeling processes like super exchange and singlet fission, 13 trivial crossing. 14 SH has been successfully applied to model charge and energy transfer, carrier recombination in various systems, but it is still limited to small systems with a rigorous scheme. Classical Path Approximation (CPA) 15 is applied with SH when studying the solid phase systems. Given its wide applications, SH is still under development to improve its efficiency and accuracy. 3 While CMD is notorious for its accuracy especially when transferring to model a new system, and ab initio MD is computationally costly, the incorporation of machine learning is promising to address these problems. Machine Learning has been widely applied in many disciplines and earns a reputation with its extraordinary performance. Machine learning is helping revolutionize the development of force field in CMD, 16 characterizing properties 17 and at the same time provides reliable and useful techniques for analysis. 18 1.1 Classical and Born-Oppenheimer Molecular Dynamics The simulation of molecular dynamics is achieved by the numerical, step-by-step integration of the classical equations of motion, and for an atomic system the motions of nuclei are described by the Lagrangian: 𝑳 = 1 2 ∑m i 𝒓 ̇ 𝒊 𝟐 𝑖 −𝑈 (𝒓 ) (1.1) 𝒇 𝒊 = 𝒅𝑳 𝒅𝒓 =− 𝒅𝑼 (𝒓 ) 𝒅𝒓 (𝟏 .𝟐 ) Here U(r) is the potential energy surface. In CMD, U(r) is an empirical formula with fitted parameters. Like for noble gas, van der Waals interaction is described by Lennard-Jones function: 𝑈 (𝑟 𝑖𝑗 )=4𝜖 [( 𝜎 𝑟 𝑖𝑗 ) 12 −( 𝜎 𝑟 𝑖𝑗 ) 6 ] (1.3) Here rij is the interatom distance, and 𝜎 and 𝜖 are length and energy parameters. And in Quantum Molecular Dynamic, U(r) is the potential energy with the ground state electronic wavefunction Ψ 𝑜 : 𝑈 =<Ψ 𝑜 |𝑯 |Ψ 𝑜 > (1.4) H is the Hamiltonian for the system. The forces are computed based on Hellmann-Feynman theorem: 4 𝒇 𝑖 =<Ψ 𝑜 | 𝜕 𝑯 𝜕 𝒓 𝒊 |Ψ 𝑜 > (1.5) Velocity Verlet19 algorithm is usually used for integration of Eq. (1.2): 𝒓 (𝑡 +∆𝑡 )=𝒓 (𝑡 )+𝒗 (𝑡 )∆𝑡 +∆𝑡 2 𝒇 𝑚 (1.6) 𝒗 (𝑡 +∆𝑡 )=𝒗 (𝑡 )+ ∆𝑡 2𝑚 [𝒇 (𝑡 +∆𝑡 )+𝒇 (𝑡 )] (1.7) Here ∆t is the timestep for the discrete integration scheme and it is usually chosen to be hundredth of typical timescales of systems, and typical values are 0.1 to several femtoseconds. This integration scheme is time-reversible and symplectic and conserves volume in phase space. Due to the Lyopunov instability,20 the integration does not generate real trajectories but close ones, and thus it is good for properties requiring statistical averaging. This basic form is in Eq. (1.1) is usually used to sample microcanonical ensembles as there are no interaction terms with environment. Extension to sampling canonical and grand canonical ensembles are achieved with so-called thermostats and barostats by introducing extended variables.1 Like with a Nose Hoover thermostat, the Lagrangian is modified as: 𝑳 = 1 2 ∑m i 𝑠 2 𝒓 ̇ 𝒊 𝟐 𝑖 −𝑈 (𝒓 )+ 1 2 𝑄 𝑠 ̇ 2 −𝑔 𝑘 𝐵 𝑇𝑙𝑛 (𝑠 ) (1.8) here s is the extended variable and Q is the effective mass term, g controls the coupling strength. With the trajectories generated by the simulations when systems reach equilibrium, various properties of the system can be calculated with theories from statistical mechanics. In linear response theory, the susceptibility can be computed from the fluctuation of related quantities, for instance, the heat capacity can be related to the fluctuation of total energy in a canonical ensemble: 𝐶 𝑣 = 1 𝑘 𝐵 𝑇 2 <∆𝐸 2 > 𝑇 (1.9) 5 And transport coefficients, which are another group of properties of great interest, like thermal conductivity, diffusion coefficients, viscosity, can be computed from their equilibrium correlation functions. For instance, the thermal conductivity 𝜆 can be computed with Green Kubo formula:21 𝜆 = 1 3𝑉 𝐾 𝐵 𝑇 2 ∫ <𝒋 (0)𝒋 (𝑡 )>𝑑𝑡 ∞ 0 (1.10) Here j(t) is the heat current which can be evaluated with basic quantities. Apart from estimating properties from equilibrate simulations, non-equilibrate molecular dynamics can be also used to compute transport properties as they usually require extremely long time simulations for evaluating like in Eq. 1.8. One way is to create systems in steady states. For example, to compute thermal conductivity, heat is added into the system and taken out from other part to create a steady heat flow, and the thermal conductivity is directly estimated from the temperature gradient. Most importantly, non-equilibrate molecular dynamics can also be used to model processes like chemical reactions,3 crystallization,22 shock,23 etc, and to study system response upon external stimulus or to reveal the mechanisms of dynamic processes. In Chapter 2, we provide a throughout study of generating layered materials with performing shock in liquid solid system with non-equilibrate molecular dynamics. An atomic view of exfoliation process is presented, and the mechanism of exfoliation process is revealed. 1.2 Non-Adiabatic Molecular Dynamics with Surface Hopping While the adiabatic approximation breaks down in systems with excited electrons, the quantum molecular dynamics described above can be modified to use excited state forces for modeling these systems. 24 However it is computationally unfeasible if the systems have dense degrees of freedom. The popular approach is still to propagate the nuclear subsystems on adiabatic potential surfaces and at the same time to allow the system to hop between different adiabatic 6 potential surfaces. The scheme is called surface hopping (SH). 10, 25-26 Here electron dynamics are described by time-dependent Schrodinger equation (TDSE): 𝑖ℏ 𝜕 Ψ(𝐫 ,𝐑 ,t) 𝜕𝑡 =𝐻 (𝐫 ,𝑹 ,𝑡 )Ψ(𝐫 ,𝐑 ,t) (1.11) where r, and R are the electron and nuclei coordinates, H is the Hamiltonian for electrons. Ψ is the wavefunction for the electrons and is usually expressed with the eigenstates of H: Ψ(𝐫 ,𝐑 ,t)=∑𝑐 𝑖 (𝑡 )Φ 𝑖 (𝐫 |𝐑 ) 𝑖 (1.12) Eq. (1.11) is rewritten as a set of equations for the expansion coefficients ci(t): 𝑖ℏ 𝜕 𝜕𝑡 𝑐 𝑖 (𝑡 )=∑(𝘀 𝑖 𝛿 𝑖𝑗 +𝑑 𝑖𝑗 )𝑐 𝑗 (𝑡 ) 𝑗 (1.13) Here 𝘀 𝑖 is the eigen energy for Φ 𝑖 and 𝑑 𝑖𝑗 is the non-adiabatic coupling (NAC) between state Φ 𝑖 and Φ 𝑗 : 𝑑 𝑖𝑗 =−𝑖ℏ<Φ 𝑖 | 𝜕 𝜕𝑡 |Φ 𝑗 >=−𝑖ℏ<Φ 𝑖 |∇ 𝐑 |Φ 𝑗 >𝐑 ̇ (1.14) When the hopping happens between specific adiabatic potential surfaces i and j, the velocities of atoms are rescaled along the non-adiabatic coupling vector <Φ 𝑖 |∇ 𝐑 |Φ 𝑗 > to conserve energy. Different SH schemes are developed for modeling the hopping processes. The basic one, fewest switch surface hopping (FSSH), gives the hopping probability for a given time duration as: 𝑃 𝑖 →𝑗 (𝑡 ,𝑑𝑡 )=∫ 2 ‖𝑐 𝑖 (𝑡 )‖ 2 Re[( 𝑖(𝑑 𝑖𝑗 ) ℏ )𝑐 𝑖 ∗ (𝑡 )𝑐 𝑗 (𝑡 )]𝑑𝑡 𝑡 +𝑑𝑡 𝑡 (1.15) And in Decoherence-Induced Surface Hopping (DISH) 11 , the hopping is sampled from Poisson process with mean waiting time as the decoherence time 𝜏 𝑖𝑗 between adiabatic states i and j. In Chapter 3, we apply the SH to study the electron dynamics in carbon nanotube and nanoribbon and reveal the origin of different relaxation timescales. In Chapter 4, we extend the 7 NAMD scheme into modeling of multi-particle processes, as the non-adiabatic coupling is for single-particle transition. And in Chapter 5, we connect unsupervised machine learning (described in the following section) with NAMD to help analyze the simulation results. 1.3 Machine Learning in Computational Material Science Two main types of machine learning (ML), unsupervised and supervised ML have been widely used in material science. They are incorporated to address problems like improving the model accuracy, reducing the computational cost, relieving the effort for developing models. One is unsupervised learning used in material science for material design and discovery. Here statistical methods are used to extract hidden patterns, information and complex dependencies from bunch of data which has similarities and differences. Properties like dielectric or viscosity have implicit dependencies on the system configurations. When we search for candidate materials with desired properties, unsupervised learning can help extract this dependency which normal theoretical methods can’t provide directly or accurately, and greatly reduce the experimental work for synthesizing and testing. For example, to search for polymers with high thermal conductivity, Wu and coworkers apply unsupervised learning assisted model and generate polymer systems verified by experiments to exhibit high thermal conductivities. 27 In Chapter 5, we apply one unsupervised learning algorithm to estimate mutual information between quantities from nuclei and electron motions, and extract critical motions of nuclei affecting the electronic dynamics. The other one is supervised learning, where the training data includes input or features, and corresponding targets or output. For example, to build a machine learning model for performing molecular dynamics, we can take system configurations as input and potential energies and/or atomic forces as output, as such we can use the integration scheme in Eq. (1.6) and (1.7). The 8 applied algorithms in material science range from simple linear regression, to Kernel ridge regression and neural network. And it has been successfully used for predicting properties like potential energy surfaces, atomization energies, forces, atomic partial charges, molecular dipoles. 28-29 Given these successful modeling and applications, machine learning is criticized for its poor interpretability. Many techniques, especially the most popular one, neural network, work like a black box, and provide few intermediate physical information. Apart from this, to have a higher accuracy, the training of these models usually asks for huge amount of data points. For example, the AN1 model is trained with tens of millions of molecular systems. 16 These data points are usually generated from first principle calculations. Generating of these data points and training these models requires lots of computation resource. Last but not least, due to the tremendous chemical compound and configuration spaces, these models usually have poor transferability and extensibility, i.e., performing poorly on systems which has different chemical bonding or elements as in training set, and on systems which are larger than the systems in training set. To improve its transparency, transferability, extensibility and at the same time to reduce the requirement of data amount, a promising method is to connect machine learning with quantum chemistry methods. In Chapter 6 and 7, we provide a detail description of semi-empirical quantum chemistry package PYSEQM and incorporate it with a specially design neural network for chemical systems. This hybrid model shows an overall better performance on predicting atomization energies and atomic forces. Furthermore, it can be extended to model excited state molecular dynamics and study a wide range of phenomena. Incorporating machine learning in computational material science is better than using traditional methods in many aspects. Machine Learning not only provides many techniques for 9 analyzing data in materials science, but also help address the problems in theoretical study. In computation material sciences, the high-level methods like Density functional theory, Coupled Clusters 30 are computationally costly and are restricted to model small systems, while the low- level methods like classical molecular dynamics are efficient but less accuracy. Connecting with machine learning can put forward this boundary and limit. On the other hand, the development of traditional methods usually requires human intuition, and thus is less likely to be adapted by other researchers in the science community. For example, the development of reactive force field requires human step-by-step guides. 31 As a comparison, the development of force fields with neural network automates these processes, and introduce less bias from human. 1.4 Thesis Organization This thesis is organized as follows, and each chapter is more or less self-contained. First, we discuss two typical examples on modeling nanoscale systems in Chapter 2-3. In Chapter 2, we employ classical molecular dynamics to study the mechanism in liquid phase exfoliation with a layered system MoS2. With the force field model developed to properly describe the interaction between solid and solvent systems, the molecular dynamics simulations allowed us to have an insight view during onset of separation of atom layers and to elucidate the role of shock applied on the solvent. In Chapter 3, we study and compare the carrier dynamics in carbon nanotube and nanoribbon via non-adiabatic molecular dynamics and show nanoscale systems can exhibit distinct properties through different structures even they have same local bonds. We come to the method development in Chapter 4-7. In Chapter 4, we present a method with time-dependent density functional theory scheme to model Auger-type processes in nanomaterials. In nanoscale system, the strong electron-electron interaction from the spatial confinement will lead to multi-particle scattering, like Auger processes, and this scattering may dominate carrier dynamics. We 10 incorporate a proper description of electron-electron interaction, combine with surface hopping, to model this Auger process. In Chapter 5, we adapt an unsupervised machine learning technique into the analysis of lattice deformation and vibration impacts on the carrier dynamics and demonstrate in a hybrid halide perovskite system. The new analysis method provides a quantitative estimation and comparison of lattice impact from the complex dependency between the lattice motions and carrier dynamics. In Chapter 6 and 7, we provide an early stage study of incorporating domain knowledge from quantum chemistry into artificial neural network. This modern machine learning technique is famous for its wide application and high accuracy and efficiency, while it is notorious for its interpretability and extensibility when applying in material science. In Chapter 6, we describe the implementation of one quantum chemistry method with machine learning framework PyTorch and benchmark its performance with modern graphic process unit architectures. In Chapter 7, we apply the developed package and link with a hierarchical neural network to enhance its accuracy and transferability on nanoscale systems. Finally, in Chapter 8, we summarize and outline future directions of modeling nanoscale systems and discuss the range of problems which these developed methods can be applied to study with proper and useful accuracy and efficiency. 11 Chapter 2: Molecular Simulation of MoS2 Exfoliation Reprinted with permission from Sci. Rep. 8, 16761 (2018), Copyright 2018, Springer Nature. Authors: Guoqing Zhou, Pankaj Rajak, Sandhya Susarla, Pulickel M Ajayan, Rajiv K Kalia, Aiichiro Nakano, Priya Vashishta Contributions: R.K.K., A.N., and P.V. designed the molecular dynamics simulation project and G.Z. carried out the simulations and analyzed the results. Visualizations were done by G.Z. and P.R. Experiments were performed by S.S. and P.M.A. 2.1 Introduction Liquid-phase exfoliation (LPE) 32-38 is a highly promising approach to large-scale production and dispersion of a wide variety of layered materials (LMs). It affords facile processing of individual nanosheets, which can be deposited on surfaces or combined into free-standing films 37 , and vertical or horizontal stacks. LPE has been used to create novel multi-ferroic materials for photoconducting cells 38-40 , p-n junctions, field-effect transistors 41-42 , and memory devices 43 using stacks of layered perovskites. Integration of LMs using LPE has potential applications in large-area electronics 44 and inkjet printing 45 . A wide variety of transition metal dichalcogenides (TMDCs), metal oxides, and perovskites have been exfoliated into 2D layers by electrochemical, sonication and shear methods 32, 37 . Here, we will focus on the sonication approach in which a bulk solid is suspended in a suitable solvent and exfoliated into atomically thin LMs by cavitation phenomenon. Since LMs are characterized by strong in-plane covalent bonds and weak out-of-plane interactions, it is possible to exfoliate LMs by weakening the interlayer van der Waals interaction using ultrasonic cavitation 33 . It is desirable to keep ultrasonic intensity low in order to avoid sonolysis and defects in 2D materials and to prevent radicals in solvents which may affect dispersion of LMs. Solvents 12 play a critical role in efficient production of 2D materials by LPE. Experimental measurements 46- 47 of interfacial energy, Hildebrand solubility and Hansen parameter are commonly used to guide the selection of solvents for exfoliation of LMs. Solvents with weak volatility (N,N- dimethylformamide and N-methyl-2-pyrrolidone) and low boiling points (propanol, chloroform) have been successfully used to exfoliate LMs 32, 46 . Despite a great deal of experimental work, there is very little understanding of atomistic mechanisms underlying LPE. The motivation for the joint experimental and simulation work reported in the paper is to unveil the atomic mechanism of liquid-phase exfoliation and thus facilitate the synthesis of atomically thin layered materials (LMs). Shock exfoliation of LMs by bubble collapse mimics experimental conditions, and experimentalists are trying to optimize the conditions for liquid-phase exfoliation by choosing suitable solvents and shock intensity so that single sheets of defect-free LMs are produced. We have performed molecular dynamics (MD) simulations 48-49 in which a MoS2 crystal is suspended in a solvent of water and isopropanol (IPA) containing a cavitation bubble. The system is subjected to a planar shock which initiates a chain of events in the solvent, culminating in the exfoliation of MoS2 into nanosheets. As the shock wave propagates through the solvent, the density of solvent increases to 1.6 g/cc and the number of nearest neighbors of a water molecule increases from 4 to 8, indicating an ice VII like motif. Water in the compressed solvent is not frozen: on the contrary, the self-diffusion coefficient of H2O molecules normal to the direction of shock propagation is increased by 60%. The shock wave impact collapses the cavitation bubble and generates a high-speed nanojet in the solvent. The nanojet impact generates large shear stresses (~10 GPa) on the MoS2 surface and the surface temperature goes up to ~3,000K. These large shear stresses and elevated temperature initiate exfoliation of MoS2, and shock waves reflected from MoS2 surfaces enhance exfoliation. We have 13 performed LPE experiment with IPA and DI water, which verify the conclusion from the simulations. 2.2 Simulation Setup Figure 2.1(a) shows an initial configuration of a simulation in which an MoS2 solid (yellow and pink spheres) is immersed in a 1:1 mixture of H2O and IPA. (For the sake of clarity, only the lower half of the MD box is shown in the figure.) Cavitation is introduced after equilibrating the system under ambient conditions. The ratio of the bubble radius to the shortest distance between the bubble center and MoS2, i.e. the stand-off parameter S, plays a critical role in exfoliation. We performed simulations for several values of S ranging between 1 and 2 and observed exfoliation in the range 1.1 < S < 2.0 with particle velocity Vp = 3.0 km/s. Here we will present results for a bubble of diameter 9.4 nm and S = 1.14. Additional results are presented in the supplementary information. 14 Figure 2.1: Initial configuration of the shock-induced exfoliation simulation and formation of ice VII motif. (a) shows the initial setup of the exfoliation simulation. Bulk MoS2 is represented by pink (Mo) and yellow (S) sheets. MoS2 is immersed in a solvent, consisting of water and IPA molecules (1:1 ratio by weight). For clarity, only 2% of the solvent molecules (Oxygen: red, Carbon: cyan, Hydrogen: white) are shown in the lower half of the MD box. The solvent contains a nanobubble of radius R = 4.7 nm. The stand-off parameter, d/R = 1.14, where d is the distance between the bubble center and the closest MoS2 surface. (b) shows the radial distribution function for oxygen-oxygen in water (red) and oxygen in water and the center of mass of IPA (blue). (c) shows one water molecule with 8 nearest neighbors. Six of them are H2O and two are IPA molecules. In all simulations, shock is generated by a momentum mirror placed normal to the z direction just outside the MD box. The solvent, MoS2 and bubble are moved towards the momentum mirror with a constant speed Vp at time t=0,. When the solvent molecules cross the mirror, their momenta in the z direction are reversed which creates a planar shock wave in the solvent propagating away from the mirror. Using this approach, we first calculated the Hugoniot 15 (shock-wave velocity Vs as a function of Vp) of the solvent without MoS2 and found that it was similar to the Hugoniot of pure water; see Figure S2.2 in the supplementary information. 2.3 Results and Discussion 2.3.1 Shock Impact on Solvent Shock simulations are performed at several particle velocities in the range of 0.5 - 4.0 km/s. Here we present results for the shock velocity Vs = 7.4 km/s corresponding to the particle velocity Vp = 3.0 km/s. A movie of exfoliation is in the supplementary material. Under these conditions, the pressure in the solvent rises to 10.5 GPa and the density of water increases from 0.95 g/cc to 1.59 g/cc. This high compression has a dramatic effect on the structure of water. Figure 2.1(b) shows the radial distribution function go-o(r) for oxygen atoms of water molecules. Here the first peak is located at 2.68 Å and the first minimum at 3.75 Å, whereas in pure water under ambient conditions the first peak and first minimum are at 2.76 Å and 3.34 Å, respectively. The average number of nearest neighbors of water molecules, calculated from the area under the first peak in Figure 2.1(b) with a cutoff distance of 3.75 Å, is 8 as opposed to 4 in pure water under ambient conditions. The O-O-O bond-angle distribution for water molecules in the high-density solvent (HDS) peaks around 56 ◦ and the O-O-IPA distribution peaks around 65 ◦ . These results indicate that the structure of water in the HDS is similar to that of ice VII in that both of them have the same density (1.6 g/cc) and number of nearest neighbors (8), see Figure 2.1(c). However, the differences in bond-angle distributions reflect disorder in the first solvation shell of water in the HDS compared to the first solvation shell of ice VII 50-51 . The MD results for the structure of water in the HDS are in good agreement with Dolan et. al.’s shock-wave experiment 51 on pure water. They observed rapid freezing of water for Vp between 0.5 km/s and 2.0 km/s under isentropic and ramp-wave compression. Freezing occurred 16 within a few nanoseconds above a critical value of pressure (7 GPa) irrespective of the peak pressure generated by shock. Their data clearly show that ultrafast homogenous nucleation of ice VII is feasible under shock compression above 7 GPA. In the MD simulation, we do not observe complete freezing of H2O into ice VII because the timescale of the applied shock (~ 20 ps) is much shorter than the timescale (a few ns) in the experiment and also because of the presence of IPA molecules. However, we do observe a significant change in the dynamics of water molecules in the HDS: the self-diffusion coefficient of water, calculated from mean-square displacements normal to the direction of shock-wave propagation, is larger (3.710 -5 cm 2 /s) than that of pure water under ambient conditions (2.410 -5 cm 2 /s). The IPA molecules in the HDS diffuse much more rapidly than under ambient conditions: the self-diffusion coefficient for the center-of-mass motion of IPAs in the 50%wt mixture is 2.310 -5 cm 2 /s in the HDS and 1.410 -6 cm 2 /s under ambient conditions. 17 2.3.2 Shocked Induced Bubble Collapse Figure 2.2: Snapshots of the collapsing cavitation bubble. (a) shows the initial configuration of the cavitation bubble and MoS2. The bubble is represented by a shell of liquid molecules on the bubble surface. The standoff parameter is 1.14. At t = 0.2 ps, the shock wave hits the proximal side of the bubble. (b) – (d) are snapshots showing changes in the shape of the collapsing bubble due to the shock wave at time t =0.65, 0.85, 1.1 ps. The shock wave is represented by the change in the liquid density (half of the density is shown in the figures). (e) shows changes in the surface area and volume as a function of time while the bubble is collapsing. The surface area and volume are normalized to their respective initial values. The bubble collapses at t = 1.5 ps. Under the impact of the shock wave, the bubble begins to shrink because the surface tension of the bubble cannot provide enough restoring force to balance the shock-wave compression. Snapshots in Figure 2.2 show a time sequence of changes in the shape and size of the nanobubble resulting from the shock impact. As more and more solvent molecules enter the bubble, the proximal side of the bubble changes from spherical to ellipsoidal. The shape of the shock front also changes during bubble shrinkage: the front loses planarity because the solvent molecules entering the bubble have different velocities than the shock-front velocity Vs. The front regains planarity after the bubble disappears. The bubble collapse time – the elapsed time between the onset of bubble shrinkage and complete bubble collapse – is 1.5 ps, see Figure 2.2(b). It is in close 18 agreement with the Rayleigh formula for bubble collapse time, 𝜏 =0.45𝐷 √ 𝜌 ∆𝑃 , (2.1) where D is the initial diameter of the bubble,  is the fluid mass density (HDS in our case), and P is the pressure difference across the bubble surface. Substituting D = 9.4 nm,  = 1.59 g/cc and P = 10 GPa in Eq. (2.1), we find  = 1.7 ps. It is remarkable that this estimate agrees so well with the MD result even though the effects of viscosity, surface tension and non-uniformity of the solvent near the bubble surface are ignored in the Rayleigh formula. 2.3.3 Onset of exfoliation At the onset of the bubble collapse, we notice sudden increases in the translational kinetic and rotational energies of solvent molecules at the shock front. These energy jumps are caused by solvent molecules entering the proximal side of the bubble. Velocity streamlines of these molecules are focused towards the bubble center in the form of a high-speed nanojet. Figure 2.3(a) shows the center-of-mass velocity streamlines of molecules in the nanojet at t = 1.75 ps after the bubble collapse. The nanojet length and width are 6.5 nm and 3 nm, respectively. At the tip of the nanojet, the velocity is 6.1 km/s and the pressure around 20 GPa, which is close to the pressure estimated from the jump condition, P – P0 =  VpVs = 21 GPa. The nanojet length increases linearly with the initial diameter of the bubble 52 . Experiments 53-54 show that this linear relationship also holds for micron-to-millimeter size bubbles. The nanojet generates a vortex (see the inset in Figure 2.3(b)) whose angular velocity, calculated from the stream velocity Ω ⃗⃗ =∇×v ⃗ , ranges between 5 and 15 ps -1 . The nanojet hits the MoS2 surface at t = 1.75 ps after the bubble collapse. The impact causes pit formation on the MoS2 surface (indicated by the red region Figure 2.3(c)), resulting in 19 an 11% volume reduction at t = 3ps. The pit is 3 nm wide and 1 nm deep. At t = 3.7ps, the convex hull volume of MoS2 expands by 20%. Figure 2.3: Shows cavitation bubble collapse giving rise to a nanojet and a vortex. (a) Snapshot taken at t = 1.75 ps shows velocity streamlines of the nanojet resulting from the bubble collapse. The color represents the magnitude of the stream velocity, which ranges between 3 and 6 km/s. The nanojet length and width are 6.5 nm and 3 nm, respectively, and the pressure at the tip of the nanojet is around 20 GPa. (b) The inset shows a vortex generated by the nanojet. The nanojet impact on MoS2 creates a 3 nm wide and 1 nm deep pit at t = 3 ps. (c) von Mises shear strain distribution in the pit region initiates exfoliation of MoS2. 20 Figure 2.4: Pressure, temperature and shear stress distributions in MoS2 after the nanojet impact. (a) At t = 2 ps, the pressure is 40 - 50 GPa in the red region; 10 - 20GPa in the green/yellow region; and around 0 GPa everywhere else. (b) shows that the pressure is around 10 - 20GPa over almost the entire MoS2 at t = 3 ps. (c) shows that the pressure is released at t = 40 ps. (d) Shows the antisymmetric distribution of shear stress Sxz at t = 2 ps, (-10 to -5 GPa in the blue region, and 5~10 GPa at the center). (e) shows Sxz spreads to initiate exfoliation at t = 3 ps. (f) shows that the shear stress is released at t = 40 ps. (g) Temperature in the pit (red) is around 3,000K, whereas the rest of the sample is around 300 K. (h) shows that the high temperature region becomes 7 nm wide and 1.5 nm deep. Here the temperature in the green and yellow regions is around 1,500 K. (i) Shows that the temperature in MoS2 is uniform and around 1,600 K after exfoliation. 21 The nanojet impact has a dramatic effect on exfoliation of MoS2. The time evolution of exfoliation is related to pressure, shear stress and temperature distributions in MoS2, see Figure 2.4. The kinetic energy imparted by the nanojet significantly increases the pressure in the pit region of the MoS2. At t = 2 ps, the instant pressure in the 4 nm wide and 1 nm deep pit (red region in Figure 2.4(a)) varies between 40 and 50 GPa. The pressure in the green and yellow regions around the pit ranges between 10 and 20 GPa, whereas the rest of the MoS2 is still at zero pressure. The pressure begins to drop as the pit region expands to 7 nm in width and 4 nm in depth at t = 2.3 ps. At 3.0 ps the pressure is between 20 and 30 GPa in almost all of MoS2 (see Figure 2.4(b)), and after t = 7 ps the pressure drops to 0 GPa. Figure 2.4(c) shows the pressure distribution after MoS2 exfoliation at t = 40 ps. Temperature fluctuations closely follow changes in the pressure in MoS2. The local temperature (calculated with peculiar velocity) in the pit ranges between 2,500 and 3,500K at time at t = 2ps. Away from the pit, the MoS2 remains at room temperature (indicated by the blue region in Figure 2.4(g)). At t = 3ps, the pit region cools off slightly due to expansion and the temperature distribution ranges between 2,500 and 3,000K. At t = 4ps, the temperature in MoS2 drops to 1,000K except at the surface impacted by the nanojet where the temperature is 1,500K. The temperature fluctuations persist until the very end of exfoliation, see Figure 2.4(i). The nanojet impact also generates large shear stresses, which initiate exfoliation of MoS2 into nanosheets. Figure 2.4(d) shows that the shear stress component Sxz in the red (compressive) and blue (tensile) regions fluctuates between 5 and 10 GPa, whereas in the rest of the MoS2 sample the shear stress is negligible. (Note, the shear-stress distribution is antisymmetric along the central plane at x = 100 Å.) A similar pattern is observed in the shear stress component Syz, which is antisymmetric along the plane y = 100Å (see Figure S2.2 in supplementary information). These 22 large shear stresses initiate exfoliation of MoS2 layers. Figure 2.4(e) shows that bulk MoS2 has partially exfoliated and the shear stress has dropped to 5 GPa at t = 3 ps. The shear stress is released after exfoliation, see Figure 2.4(f). 2.3.4 Exfoliation with repeated impact of shock wave The exfoliation is significantly enhanced every time the shockwave propagates through MoS2; see Figure S2.3 in the supplementary information. The instantaneous temperature of MoS2 increases to 1,300K during the first encounter with the shock wave between 1.75 and 3 ps. Subsequently, the convex hull volume increases by 20% as 7,500 H2O and 1,700 IPA molecules flow inside the galleries of MoS2. The temperature of the sample decreases to 1,100 K at t = 12 ps and about 73% of the solvent remains in the MoS2 layers. The shock wave is reflected from the back end of the MD box (z = 0) at t = 8 ps and when the release wave hits MoS2 at t = 12 ps, the temperature of MoS2 increases to 1,400 K and the solvent content in the galleries of MoS2 increases with the addition of 3,000 water and 800 IPA molecules. Eighty percent of the solvent remains between the MoS2 layers after the passage of the release wave. As the shock wave reaches the opposite end of the box (z = 28.7nm) and reflected again at time t = 18 ps, the temperature of MoS2 increases to 1,650 K at t = 23 ps, and an additional 5,700 H2O and 1,500 IPA molecules enter the galleries of MoS2. Figure S2.3 in supplementary information shows solvent molecules between MoS2 nanosheets at t = 40 ps. Only 10% of the molecules chosen randomly are shown here. There are 13,500 water and 3,300 IPA molecules inside the sample at t = 40 ps. The volume of the convex hull continues to increase as more nanosheets exfoliate, see Figure S2.3(b). The swelling parameter of MoS2 plateaus at 2.0. The number of solvent molecules in the galleries of the exfoliated nanosheets increases rapidly when shock waves hit the MoS2 at t = 1.75, 12, and 18 ps; see Figure 23 S2.3(c). Some of these solvent molecules diffuse out of the galleries. Overall, there is a dynamic balance between the number of solvent molecules entering and leaving the galleries of MoS2 sheets. Figure S2.3(c) in the supplementary information shows that the number of water molecules flowing in and out of the MoS2 sheets is larger than the number of IPA molecules because MoS2 is hydrophobic and water diffuses more rapidly than IPA. 2.3.5 Experimental Results Figure 2.5: Exfoliated MoS2 with solvent between the nanosheets. (a) Commercial MoS2 in IPA+DI water solution (i) before and (ii) after exfoliation. (b) Raman spectrum of the exfoliated MoS2 flakes. (c) Snapshots of MoS2 flakes being exfoliated in IPA+ DI water mixture. A small bubble formation takes place at the interface of MoS2 and solvent. The bubble expands after 90s and collapses at 120s, taking MoS2 powder along and dispersing it in the solvent. The color change in the liquid is an evidence of this phenomenon. To verify the theoretical claims, experimental exfoliation of MoS2 was performed in IPA and DI water mixture (see Methods section). Figure 2.5(a) shows the black colored MoS2 powder in IPA+DI Water solvent before exfoliation. After exfoliation, the color of the solvent changes from transparent to light green due to the presence of 2D MoS2 flakes dispersed in the solvent. A 24 slow-motion video was captured during exfoliation to examine different stages of exfoliation. Figure 2.5(c) shows snapshots of MoS2 powder at different stages of exfoliation. First, a bubble forms at the interface of MoS2 and solvent as a result of sonic wave propagation in the solvent. The bubble expands after 30s and collapses after 60s, resulting in shock wave inside the solvent and MoS2. The bubble also carries MoS2 with it, which is indicated by the local change in color of the solvent when the bubble bursts. The process repeats itself several times, resulting in uniform dispersion of 2D MoS2 sheets in the IPA+DI water solvent. 2.4 Conclusions In conclusion, MD simulations reveal atomistic processes underlying liquid-phase exfoliation (LPE) of LMs by sonication. Cavitation phenomenon underlies LPE experiments in which sonication probes are used to generate cavitation bubbles. We find that bubble collapse is a highly energetic process, and this is corroborated by experimental evidence of hot-spot formation in solvents 55 . Experiments 56 indicate that temperatures in hot spots can be as high as 5,000 K. Our simulation shows that shock-induced bubble collapse results in the formation of a high-speed jet whose impact on an MoS2 surface initiates the exfoliation process. The nanojet impact raises the MoS2 surface temperature to 3,000 K and pressure to 20 GPa, and exerts a shear stress of 10 GPa. Exfoliation of MoS2 is initiated by this combination of large shear stress and temperature increase and, finally, exfoliation is significantly enhanced by repeated interactions of MoS2 with release waves resulting from the reflection of shock waves in the system. The experiment has helped us validate the simulation. Experiments reveal that sonication causes cavitation and the collapse of cavitation bubbles generates shock waves in the solvent. Bubble creation and collapse are captured with high resolution camera and Raman spectra show that bulk MoS2 is exfoliated and flakes of MoS2 are observed. Our simulations are in accord with 25 experimental observations. The simulation results - stress and temperature distributions, size and placement of bubbles, nature of solvent and solvent concentration - can help experimentalists with the optimization and scaling of exfoliation yield. 2.5 Method In the MD simulation, an MoS2 crystal of dimensions (9.8 nm) 3 is suspended in an H2O/IPA mixture (1:1 ratio by weight) containing a spherical cavitation bubble of radius 4.7 nm. The system contains 10 6 atoms and its dimensions are 19.7nm×19.7nm×28.7nm in the x, y, and z directions, respectively. Periodic boundary conditions are applied along x and y and fixed boundary condition in the z direction. The simulation is carried out with a combination of force fields: TIP4P/2005 57 for water, REBO potential 58-59 for MoS 2 , and OPLS-AA force field 60 for IPA. (The OPLS-AA is commonly used for organic molecules.) The interaction between MoS2 and water is described by a combination of Lennard-Jones (LJ) and electrostatic potentials with force-field parameters taken from Luan et al. 61 , and we apply the Lorentz-Berthelot combination rule to parameterize LJ interactions between H2O, MoS2 and IPA molecules. Electrostatic forces and energy are calculated with the PPPM method, and slab correction 62 is applied to allow for fixed boundaries in the z direction. Force fields are validated by experimental data on contact angles and surface tensions 46, 63 , see the supplementary information Section 1. We use the Velocity-Verlet integrator with constraints on O-H bond and H-O-H bond angles of water molecules imposed with the SHAKE algorithm 19 . All simulations are done with open source package LAMMPS 49 and the visualization is done with OVITO and ParaView 64-65 . The system is relaxed for 250 ps at 300 K with a time step 0.5 fs and then a spherical nanobubble of radius 4.7 nm is created in the solvent. Exfoliation of MoS2 depends critically on the stand-off parameter, i.e., the ratio of the distance between the bubble center and the nearest 26 MoS2 surface to the bubble diameter. The stand-off parameter was varied between 1 and 2 to determine the optimum value for exfoliation. The system was subjected to planar shock in the z direction using a momentum mirror. The particle velocity was varied between 0.5 and 4.0 km/s and the time step was reduced to 0.1 fs during shock. To calculate thermal and mechanical properties in the system, the MD box is divided into 20 × 20 × 58 voxels and first the center-of-mass velocity of each voxel, 𝑣 𝑗 =∑ 𝑚 𝑘 𝑢⃗ 𝑘 /∑ 𝑚 𝑘 is computed and subtracted from the velocity of each atom inside the voxel to get the thermal velocity 𝑣 𝑘 ,𝑗 . The instantaneous “temperature” distribution in the system is determined by calculating the kinetic energy of each voxel: 1 2 𝑁 𝑓 ,𝑗 𝑘 𝐵 𝑇 𝑗 =∑ 1 2 𝑘 𝑚 𝑘 𝑣 𝑘 ,𝑗 2 , (2.2) where kB is the Boltzmann’s constant, 𝑚 𝑘 is the mass and 𝑣 𝑘 ,𝑗 is the velocity of the k th atom in the j th voxel, and 𝑁 𝑓 ,𝑗 is the number of degrees of freedom in the j th voxel. The summation is over all the atoms in the voxel.   The pressure distribution in the system is calculated from the virial stress tensor for each voxel j: 𝑆 𝑗 𝛼𝛽 = 1 𝑉 ∑[−𝑚 𝑖 𝑣 𝑖𝑗 𝛼 𝑣 𝑖𝑗 𝛽 + 1 2 ∑(𝑟 𝑘𝑗 𝛼 −𝑟 𝑖𝑗 𝛼 )𝐹 𝑖𝑘 𝛽 𝑘 ] 𝑖 ∈𝑉 , (2.3) where V is the volume of a voxel, 𝑟 𝑖𝑗 𝛼 and 𝑣 𝑖𝑗 𝛼 are the Cartesian components of the position and velocity of the i th atom in the j th voxel, respectively, and 𝐹 𝑖𝑘 𝛽 is the force on atom i due to atom k. The outer summation is over the atoms in a voxel, and the inner summation is over the atoms in 27 the neighbor lists of atom i. The pressure in a voxel is given by, 𝑃 𝑗 =− 1 3 𝑇𝑟 (𝑆 𝑗 𝛼𝛽 ) . (2.4) We have also estimated the exfoliation yield by computing the accessible surface area 64, 66 of the MoS2 sample. The surface area is calculated with a sphere of radius 4.0 Å. To compute the swelling of the sample, a 3-dimensional convex hull 67 is constructed with the largest clusters of MoS2 and the number of IPA and water molecules flowing in MoS2 galleries are counted. To measure the impact on bulk MoS2, we calculate von Mises local shear strain 68 𝘂 𝑘 Mises at each atom k. We choose the initial relaxed configuration of bulk MoS2 and the current configuration to get the local transformation matrix by minimizing ∑|𝑟 𝑗𝑘 0 𝑱 𝑘 −𝑟 𝑗𝑘 1 | 2 →𝑱 𝑘 =(∑𝑟 𝑗𝑘 0𝑇 𝑟 𝑗𝑘 0 𝑗 ) −1 (∑𝑟 𝑗𝑘 0𝑇 𝑟 𝑗𝑘 1 𝑗 ) 𝑗 . (2.5) Here the summation is over the nearest neighbors of atom k, 𝑟 𝑗𝑘 0,1 is the separation of atom j and k at the initial and current configurations. The shear strain of atom k is then computed as 𝘂 𝑘 Mises = √ 𝘂 𝑥𝑦 2 +𝘂 𝑦𝑧 2 +𝘂 𝑧𝑥 2 + (𝘂 𝑥𝑥 −𝘂 𝑦𝑦 ) 2 +(𝘂 𝑦𝑦 −𝘂 𝑧𝑧 ) 2 +(𝘂 𝑧𝑧 −𝘂 𝑥𝑥 ) 2 6 , (2.6) where 𝘂 𝑎𝑏 ( 𝑎 ,𝑏 =𝑥 ,𝑦 ,𝑧 ) are the six components of the local Lagrangian strain matrix 𝛈 𝑘 of atom k, and 𝛈 𝑘 = 1 2 (𝐉 k 𝐉 k T −𝐈 ). In the experiment, isopropanol (IPA) (99.99%; Sigma Aldrich) and de-ionised (DI) water were mixed in equal proportions. 2 mg of MoS2 powder (99.99%; Sigma Aldrich) was immersed 28 in 50 ml of solvent mixture (IPA +DI water). MoS2 powder was allowed to sonicate for 48hrs in the solvent. The temperature of the bath was maintained by changing the water bath every hour. Different stages of the exfoliation process were captured by a slow-motion video at 240 fps using a telephoto lens. 29 Chapter 3: A Comparative Study of Electron –Phonon Scattering in Carbon Nanotubes and Graphene Nanoribbons Reprinted with permission from J. Phys. Chem. Lett. 2019, 10, 22, 7179-7187, Copyright 2019, American Chemical Society. Authors: Guoqing Zhou, Chao Cen, Shuyi Wang, Mingsen Deng, Oleg V. Prezhdo Contributions: O.V.P. and M.D. designed the simulation project, G.Z. performed the simulations and analyzed the results. All of them contributed to the discussion and preparation of manuscript. 3.1 Introduction Graphene has been attracting attention for over a decade because of its unique structure, and extraordinary electronic, photonic, electric, mechanical and chemical properties. 69-72 With the thickness reduced to the single atom limit, and the inversion and time-reversal symmetries, the two-dimensional system exhibits a Dirac cone in its electronic structure, which leads to a vanishing bandgap and massless charge carriers with very high mobilities. 69 However, the gapless band structure hinders graphene applications in optics, electronic and photovoltaics, which require semiconducting properties. In order to create a semiconductor bandgap, one can apply the idea of quantum confinements, giving rise to two types of one-dimensional graphene-based nanostructures, carbon nanotubes (CNT) and graphene nanoribbons (GNR). CNTs were characterized 73 a decade before the isolation of single layer graphene. 69 Conceptually, a CNT can be viewed as a graphene derivative in which the atom-thick sheet is rolled up and connected along a particular lattice vectors. 74-75 Depending on the lattice vector, CNTs can be metallic and semiconducting, with varying diameters, energy gaps, stiffness and other properties, forming a rich family of structures. 30 Synthetically, CNTs are created by various methods, including arc discharge, lasers and catalyzed chemical vapor deposition. 74, 76-77 These tubular structures can form bundles or can be isolated into multi-walled and single-walled CNTs (SWCNTs). Current synthesis and purification technologies have a high degree of selectivity for CNT size, type and chirality. 78 Planar shaped GNRs can be viewed as nanometer-wide strips of graphene, or as conjugated polymers made of wide repeat units. GNRs can be fabricated top-down by lithographic etching of graphene and CNT unzipping, as well as bottom-up through on-surface or solution-based synthesis. 79-80 By tailoring the chirality and diameters of CNTs, as well as widths and edges of GNRs, one can control bandgaps and design many fruitful structures, 81-82 leading to a wide range of applications in electronics, 83-84 electrochemical energy storage, 81, 85-86 hydrogen storage, 87-89 photovoltaics, 90-91 electrochemical sensing, 92 etc. Applications of CNTs and GNRs depend strongly on their ability to transport charges. In turn, charge transport is governed by charge scattering processes, among which electron-phonon scattering constitutes the main mechanism of charge and energy losses to heat. Multiple time- resolved spectroscopic experiments and a few simulations have been carried out in order to characterize charge scattering in CNTs and GNRs, including both carrier-carrier and carrier- phonon scattering. 93-109 It has been demonstrated that energetic charge carriers in CNTs couple to the high frequency G-modes, while charges with lower energies couple of lower frequency radial breathing (RBM) vibrations. Atomic and boundary defects, as well as CNT ends are known to contribute significantly to nonradiative charge-phonon relaxation, since CNTs are very good conductors and charges can easily reach defect sites. Charge transport in GNRs is strongly influenced by edges, since edges are susceptible to distortion and facilitation formation of 31 localized charge carriers, known as polarons, which exhibit different transport properties than free carriers. 110-111 Particularly interesting is the work of Bonn and co-workers, 93 since it provides a direct comparison between a CNT and a GNR of similar size, using the same experimental tools. Such comparison gives an opportunity to understand and characterize the similarities and differences between these seeminly analogous nanostructures, which can be converted between each other. In particular, Bonn and co-workers measure time-resolved photoconductivity of CNTs and GNRs, and show that electron-phonon scattering events are more frequent in GNRs, and as a result, the conductivity decays several times faster in GNRs than CNTs. These, 93 as well as many other 94-99 time-resolved experiments are most closely mimicked by nonadiabatic (NA) molecular dynamics (MD) techniques, developed initially 112 in the 1960s and 1970s to study gas phase and surface scattering processes, and adapted later 10, 113-115 to condensed phase systems. Our group developed several NAMD approaches 11, 13, 26, 115-117 designed specifically for nanoscale systems, and implemented them within time-dependent density functional theory (TDDFT). 15, 118-119 In this letter, we employ large scale ab initio NAMD to demonstrate that charge-phonon scattering is much faster in a GNR than a CNT of similar size and bandgap. Even though both systems are formally quantum confined derivaties of graphene and exhibit similar electronic structure at 0 K, their properties differ significantly at room temperature. GNRs are much more flexible than CNTs and undergo significant geometry distortions at ambient conditions. The distortions localize electron and hole wavefunctions, and enhance both inelastic and elastic charge- phonon scattering, which determine charge momentum, energy and lifetime. The large inharmonicities exhibited by GNRs, as well as the presence of light H atoms at the edges, broaden the spectrum of vibrational motions that couple to charges. The NA electron-phonon coupling is 32 an order of magnitude stronger in the GNR, and the nonradiative electron-hole recombination is significantly faster. Nearly identical at the first glance, CNTs exhibit features associated with bulk materials, while GNRs show properties that are much more molecular. The reported results show good agreement with the time-resolved photoconductivity and photoluminescence measurements, provide a detailed atomistic picture and unique insights into the charge-phonon dynamics in graphene nanostructures, and suggest that controlling the rigidity of one-dimentional nanostructures can be used to tune their performance. 3.2 Computational Details Figure 3.1. Simulation cells showing optimized structures of (a) (6,5) single wall carbon nanotube (SWCNT), and (b) graphene nanoribbon (GNR). The SWCNT has 364 carbon atoms, and the GNR has 352 carbon and 132 hydrogen atoms. Electronic density of states for (c) SWCNT and (d) GNR. The study focuses on the two graphene nanostructures investigated experimentally by Jensen et al., 93 namely, the (6,5) SWCNT and GNR, as shown in Figure 3.1a,b. For computational 33 efficiency, the -C12H25 alkyl side chains in the GNR are replaced with the short methyl groups, - CH3. The (6,5) SWCNT is represented by a 111 unit cell containing 364 C atoms. The simulation cell of the GNR contains 11 unit cells, 1111, with the total of 352 C atoms and 132 H atoms. The system sizes are chosen to be similar to facilitate direct comparison of the simulation results, such as densities of states (DOS) and phonon modes involved in the nonradiative relaxation. 10 Å of vacuum are added in the two directions perpendicular to the longitudinal axis of the SWCNT and GNR, in order to eliminate spurious interactions between periodic images, and to represent isolated systems. The geometric structure optimization and adiabatic MD simulations are performed with the Vienna Ab initio Simulation Package (VASP), 120 employing the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional 121 and the projector augmented wave (PAW) method. 122 The - point is used to sample the Brillouin zone of both systems, since their atomistic representation is already large, involving over 350 C atoms in each case. The optimized structures of the SWCNT and GNR have the bandgaps of 0.97 eV and 1.51 eV, respectively, relatively close to the experimental values of 1.17 eV and 1.88 eV. 93 The bandgaps are underestimated due to the well- known self-interaction error of pure DFT functionals, such as PBE. Subsequently during the NAMD simulations we scale the PBE bandgaps to the experimental values. After the geometry optimization, the systems are heated to 300 K with the Nose-Hoover thermostat, and thermalized for 2 ps. Then, 5 ps microcanonical trajectories are generated with a 1 fs timestep. The trajectories are used to perform NAMD using the in-house version of the PYthon eXtension for Ab Initio Dynamics (PYXAID) code, 15, 119 interfaced with VASP. 120 The intraband charge-phonon relaxation is studied using Fewest Switches Surface Hopping (FSSH), 10 which is the most popular NAMD approach. The nonradiative electron-hole recombination is investigated 34 using Decoherence Induced Surface Hopping (DISH), 11 which incorporates quantum decoherence effects and produces hops at decoherence events using the standard quantum mechanical probabilities, rather than ad hoc transition probabilities such as in FSSH. FSSH is suitable to model quantum dynamics in which transitions either occur locally in time or space, e.g. near avoided crossings, or are faster than loss of quantum coherence, such as during intraband relaxation involving dense manifolds of states. Inclusion of decoherence effects via DISH, or other methods such as decoherence corrected FSSH, 101, 113 is required to model transitions that are much slower than decoherence time, 123 e.g. transitions across large energy gaps leading to electron-hole recombination. The classical path approximation (CPA) is implemented for both FSSH 15 and DISH, 119 under the assumption that thermal fluctuations in the atomic coordinates are larger than geometry changes due to photoexcitation. The CPA is essential for modeling systems involving hundreds of atoms over picosecond timescales. A detailed description of the simulation methodology can be found in references, 15, 119 which are based on the earlier work. 118, 124 The approach has been used successfully over a broad range of systems, 124-136 showing good agreement with corresponding experiments and providing a detailed atomistic description of complex quantum dynamics in nanoscale materials. 3.3 Results and Discussion The geometric and electronic structure of the optimized SWCNT and GNR exhibit many similarities, Figure 3.1. Both systems are periodic in one-dimension. Their bandgaps are on the order of 1 eV, larger for the GNR. The valence and conduction band DOS show analogous structure in the two systems, increasing away from the bandgap and showing peaks, known as van Hove singularities in SWCNTs. The main qualitative difference between the electronic states of the SWCNT and the GNR resides in the presence of light H atoms in the GNR. However, these 35 atoms do not contribute to the electronic states within the relevant energy range, since both SWCNT and GNR are -conjugated systems. The properties obtained for the optimized geometries suggest that the SWCNT and the GNR should behave very similarly, and that the charge carrier lifetime should, perhaps, be longer in the GNR, because it has a larger bandgap. 0 K 300 K SWCNT Structure HOMO LUMO GNR Structure HOMO LUMO Figure 3.2. Geometric structure, and HOMO and LUMO density distributions of the SWCNT and GNR at 0 K and 300 K. The perfect periodic structures of the SWCNT and GNR are distorted at the finite temperature. The distortion is much more significant for the GNR, because it is less rigid. Importantly, the HOMO and LUMO densities become notably localized in the GNR at the ambient temperature, while they remain delocalized in the SWCNT. 36 Consideration of the SWCNT and GNR properties at the ambient temperature already demonstrates significant qualitative differences, Figure 3.2. Both systems experience geometric distortions due to thermal fluctuations, however, the GNR distorts much more than the SWCNT. The GNR planarity is perturbed significantly. The SWCNT is no longer a perfect cylinder either. However, it maintains its shape. Most important for electron-phonon scattering and electrical conductivity, the GNR wavefunctions become localized, while the SWCNT wavefunctions are delocalized over the simulation cell even at room temperature. In this regard, the GNR behavior is similar to that of conjugated polymers, 137-138 which exhibit finite effective conjugation lengths at realistic conditions, while the SWCNT is more similar to traditional bulk semiconductors, which maintain periodicity and delocalized wavefunction over long distances. It is important to note that on length-scales longer than the simulation cell, SWCNTs will also experience distortions and bending, however, the curvature of such distortions and the extent of resulting wavefunction localization 137 are much larger than for the GNR. The wavefunction localization in the GNR at 300 K on the length-scale shorter than the simulation cell provides an indication on significant electron-vibrational interactions and has a strong effect on charge transport. The results of Figure 3.2 demonstrate the importance of finite temperature studies, since the observed effects are absent at 0 K, while they have a strong influence on the properties under investigation. The harmonic approximation near the equilibrium geometry, often used to model electron-phonon interactions, is also expected to fail for the GNR, since the harmonic approximation cannot capture the large-scale anharmonic motions illustrated in Figure 3.2. Electron LUMO+1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 SWCNT 193 866 752 627 557 406 367 317 265 242 207 205 206 GNR 310 227 149 135 131 125 124 121 116 114 111 170 166 37 Hole HOMO-1 -2 -3 -4 -5 -6 -7 -8 -9 -10 -11 SWCNT 116 1051 1023 824 727 437 420 377 355 328 309 GNR 232 180 171 144 135 156 148 128 121 199 194 Table 3.1. Intraband relaxation times [fs] for hot electrons and holes in the SWCNT and GNR, corresponding to the plots shown in Figure 3.3. Figure 3.3. Intra-band relaxation of hot electrons and holes to the band edges in the (6, 5) SWCNT and GNR with different initial excitations. The zero of energy is set to the corresponding band edge. The relaxation is significantly faster in the GNR than SWCNT. Next, we consider intraband charge relaxation, representing experimental photoexcitations at energies larger than the bandgap. 93 We consider 13 initial states for the electron and 11 initial states for the hole, Figure 3.3. These initial conditions represent photoexcitations up to 0.8 eV away from the band edges of the SWCNT and 0.5 eV for the GNR. The decay curves shown in 38 Figure 3.3 are fitted with exponents, 𝐴 exp(−𝑡 /𝜏 ) , in order to obtain the energy decay times reported in Table 3.1. Note that the fitting function contains a constant in front of the exponent. The constant accounts for the fact that early time dynamics is always Gaussian, or cosine as in the Rabi oscillation, giving rise to the quantum Zeno effect. 139-140 The evolution becomes exponential only when the wavefunction spreads over many quantum states, for instance, as required during derivation of Fermi’s golden rule. The simulation results show that the electron-phonon scattering is significantly slower in the SWCNT than the GNR. Generally, the intraband relaxation is faster for higher initial energies, because the DOS of both SWCNT and GNR increases away from the bandgap. The femtosecond intraband relaxation times obtained in the simulations agree with the carrier scattering times reported in. 93 Note that there is no direct correspondence between the two sets of data, since the experiments report carrier momentum scattering time obtained from photoconductivity measurements, while calculations obtain energy relaxation times and do not study conductivity, which requires a different type of simulation, e.g. NAMD coupled with the Landauer formalism or with nonequilibrium Green’s function theory. 141 The calculations show that charge carriers decay quickly to the band edges in both systems, indicating that after the initial sub-picosecond relaxation process, charges are conducted by states near the band edges. SWCNT Recombination [ps] 54.0 Dephasing [fs] 76.9 NA coupling [meV] 3.7  2.4 GNR Recombination [ps] 1.02 Dephasing [fs] 54.8 NA coupling [meV] 122  105 Table 3.2. Electron-hole recombination times, pure-dephasing times, and average absolute NA coupling in the SWCNT and GNR. 39 In addition to the intraband charge relaxation, we study the non-radiative electron-hole recombination across the bandgap. The recombination processes eliminate the charges and makes electrical conductivity decay to zero. The recombination timescales reported in Table 3.2 are obtained by fitting decay of the excited state population to exponent. Table 3.2 also reports the average absolute NA coupling, and the pure-dephasing times. Note that the uncertainty in the NA coupling characterizes its fluctuation along the MD trajectories, rather than an error. The uncertainty is on the order of the NA coupling value itself, because the NA couplings fluctuate significantly. The pure-dephasing times characterize elastic charge-phonon scattering taking place during the nonradiative charge recombination. The pure-dephasing times may be loosely compared with the charge momentum scattering times reported in ref 93 , since both timescales represent changes in wavefunction phase. Both pure-dephasing and momentum scattering are faster in the GNR. The nonradiative electron-hole recombination is significantly shorter in the SWCNT compared do the GNR, in agreement with the photoconductivity decay measurements. 93 The decay is faster in the GNR, because of the much larger NA electron-phonon coupling. The 1 ps electron- hole recombination time computed for the GNR matches well the 0.6 ps experimental photoconductivity lifetime obtained using optical-pump THz-probe spectroscopy. 93 The 54 ps recombination time computed for the SWCNT is notably slower than the measured 1.7 ps photoconductivity lifetime. However, as pointed out by Jensen et al., 93 their conductivity decay time is considerably lower than the photoluminescence lifetimes. 95 Our result is in good agreement with the latter data. It is important to note that the current simulation considers a perfect SWCNT. The excited state lifetime can decrease significantly in the presence of defects, 96, 101 such as point 40 defects, kinks and CNT ends, even when defect concentration is low, since charges can travel over long distances in CNTs, rapidly finding defect sites. Figure 3.4. Nonadiabatic coupling (NAC) between Kohn-Sham orbitals in (a) SWCNT and (b) GNR. The orbital indices -11, -6, 0 denote HOMO-11, HOMO-6 and HOMO, respectively, and the indices 1, 7, 14 are for LUMO, LUMO+6, LUMO+13. The size and the color of each square indicate the magnitude of the NAC for the corresponding orbital pair Further insights into the mechanism of charge-phonon relaxation in the (6,5) SWCNT and GNR are provided by Figures 3.4 and 3.5, reporting NA couplings and phonon influence spectra, respectively. Comparing the NA couplings, we observe that they are significantly larger for the GNR than the SWCNT, rationalizing the faster dynamics for both intraband and interband processes. Notably, with some exceptions, the NA couplings are significantly larger between nearest neighbor states, as indicated by the larger and brighter squares along the diagonals in Figure 3.4. This observation applies both to the SWCNT and the GNR. In order to rationalize the differences between the NA coupling values for the two systems, consider the expression used to compute these matrix elements 41 NAC=−𝑖ℏ<𝜙 𝑖 | 𝑑 𝑑𝑡 |𝜙 𝑗 >=−𝑖ℏ<𝜙 𝑖 |𝛁 𝑹 |𝜙 𝑗 >∙𝑹 ̇ (3.1) The above equation shows that the NA coupling is large when the velocities, 𝑹 ̇ , of the atoms supporting the electron and hole wavefunctions are large, and/or when the matrix element, < 𝜙 𝑖 |𝛁 𝑹 |𝜙 𝑗 >, is also large. The matrix element, <𝜙 𝑖 |𝛁 𝑹 |𝜙 𝑗 >, relies on overlap of the two wavefunctions. Also, it contains the nuclear gradient operator, 𝛁 𝑹 , and therefore, it characterizes how much electronic wavefunctions change due to atomic motions. The GNR is floppier than the SWCNT, undergoing large scale fluctuations, Figure 3.2. In addition, the GNR contains light H atoms. Therefore, the atomic velocities can be higher in the GNR. Even though the charges are supported by the -electron system of the C atoms, H atoms introduce high frequency vibrations, and high velocities, which may contribute to the electron- phonon couping in the GNR, see Figure 3.5 below. Considering the nuclear gradient matrix element, <𝜙 𝑖 |𝛁 𝑹 |𝜙 𝑗 >, note that GNR motions localize the electron and hole wavefunctions in the same place, Figure 3.2. The place where the wavefunction are localized changes along the MD trajectory, indicating that the wavefunctions are quite sensitive to atomic motions, especialy compared to the SWCNT. Therefore, the matrix elements, <𝜙 𝑖 |𝛁 𝑹 |𝜙 𝑗 >, are much larger in the GNR than the SWCNT. Thus, contributions of both the nuclear gradient matrix element and the atomic velocity to the NA coupling, Eq. (3.1), are larger for the GNR, rationalizing the faster charge relaxation and recombination. 42 Figure 3.5. Electron-phonon influence spectra for electron-hole recombination obtained as Fourier transforms (FT) of the bandgaps for (a) SWCNT and (b) GNR. Charges couple to more phonons in the GNR compared to the SWCNT, including higher frequency motions, because GNR is more anharmonic and contains light H atoms. Phonon densities of states for (c) SWCNT and (d) GNR computed as FTs of velocity autocorrelation functions. The GNR exhibits higher frequencies than the SWCNT due to presence of H atoms. Figure 3.5 a,b reports the electron-phonon influence spectra, also known as spectral densities, for the electron-hole recombination in the SWCNT and GNR. The spectra are computed as Fourier transforms (FT) of the bandgap fluctuations along the MD trajectories. These spectral densities characterize the phonon modes that promote the nonradiative charge recombination. For 43 reference, Figure 3.5 c,d present FTs of velocity autocorrelation functions averaged over all atoms of the same type in the two systems. The velocity FTs show all phonon modes available in the systems, including those that are not coupled to the charge dynamics. The nonradiative electron-hole recombination in the SWCNT is promoted most strongly by the G-mode around 1600 cm -1 arising from C-C stretching, as well as by the RBMs in the 300- 600 cm -1 frequency range, Figure 3.5a. Similar frequencies appear in the GNR influence spectrum, Figure 3.5b. Note that GNRs do not contain RBMs. An equivalent motion corresponds to GNR distortions, see Figure 3.2. Comparison of Figures 3.5a and 3.5b shows that the nonradiative charge recombination in the GNR is facilitated by a broader range of vibrational motions, including significant signals near 3000 cm -1 , corresponding to motions of the light H atoms. Interstingly, SWCNT RBMs contribute much more strongly to the velocity FT, Figure 3.5c, than the influence spectrum, Figure 3.5a. Even though the number of available higher frequency modes is smaller than the number of RBMs, the higher frequency modes correspond to lighter effective masses, have larger velocities, 𝑹 ̇ , and create stronger electron-phonon coupling, Eq. (3.1). The velocity FT for the C atoms in the GNR is much more even across the 0-1600 cm -1 frequency range, Figure 3.5d, compared to the SWCNT, Figure 3.5c. All these modes contribute to the nonradiative charge recombination in the GNR, Figure 3.5b. 3.4 Conclusions In summary, we have used the state-of-the-art time-domain simulation methodology developed in our group, combining NAMD with real-time TDDFT, in order to study charge- phonon scattering in the SWCNT and GNR of similar size. We have considered both intraband relaxation of hot electrons and holes, and interband electron-hole recombination, and have shown that these processes are significantly faster in the GNR than SWCNT. The NA electron-phonon coupling is stronger in the GNR by an order of magnitude because of the different rigidities of 44 these, seemingly nearly identical, systems. The electronic properties of the GNR and SWCNT are indeed very similar at 0 K, as suggested by the quantum confinement arguments. However, the GNR undergoes much more significant geometry distortions at room temperature than the SWCNT. As a result, the electron and hole wavefunctions become localized, symmetry selection rules for the electron-phonon interactions are relaxed, and the electron-phonon coupling becomes much stronger. The presence of light H-atoms at GNR edges enhances electron-phonon interactions further. Even though electrons and holes are supported by the -electron system of C atoms, H-atoms couple to carbons and contribute high frequency modes. Showing good agreement with the time-resolved photoluminescence and photoconductivity measurements, the reported results demonstrate that medium size GNRs with relevant bandgaps behave more as molecules than graphene, while SWCNTs do exhibit extended bulk-like features. The unique time-domain atomistic insights into the complex, non-equilibrium charge-phonon dynamics generated by the simulations modify the traditional view on these graphene nanostructures and suggest that their performance can be tuned by controlling not only size and chemical composition, but also structure and rigidity. 45 Chapter 4: Modeling Auger Processes with Non-Adiabatic Molecular Dynamics Reprinted with permission from Nano Lett. 2020, doi: 10.1021/acs.nanolett.0c04442, Copyright 2020, American Chemical Society. Authors: Guoqing Zhou, Gang Lu, Oleg V. Prezhdo Contributions: G.Z. and O.V.P. designed the project, G.Z. and G.L. implemented the model, G.Z. performed the simulations and analyzed the results. 4.1 Introduction Nanomaterials exhibit extraordinary thermal, mechanical and electronic properties, and find wide applications in chemistry, engineering, and materials science. Carrier confinement is among the key factors underlying great performance of nanomaterials. It is achieved by reduction of dimensionality, as in quantum dots (QD) and wires, and various layered systems, or by creation of superlattice structures. Finding applications in opto-electronics, quantum confinement reduces electron-phonon scattering and enhances charge carrier lifetimes. At the same time, strong confinement augments carrier-carrier scattering, amplifies dissipative Auger-type processes, and shortens carrier lifetimes. 6, 142-147 In an Auger process, one electron falls from a high energy level into a lower level, and the released energy excites another carrier, possibly leading to its emission. Auger processes depend on carrier-carrier interactions and can break the “phonon bottleneck”, arising in nanomaterials due to reduced electron-phonon scattering. 139, 148-149 Auger-mediated carrier losses can accelerate carrier trapping and recombination, undermining the electronic properties of nanomaterials. Various methods have been developed to inhibit the Auger processes in nano-systems, such as separation of electrons and holes through doping and introducing traps, 149-151 or modification of 46 nanostructure size and shape to reduce Auger efficiency. 144 At the same time, there are few theoretical works on modeling of the quantum dynamics of Auger processes, because of the expensive and time-consuming computation of electron-electron interactions. Most theoretical works compute Auger transition rates with Fermi’s golden rule based on the electron-electron interactions from static configurations. Such perturbative approach is computationally efficient, given the dense manifold of states and a tremendous number of Coulomb matrix elements. 144-145, 152-156 At the same time, Coulomb interactions can be strong under quantum confinement, challenging the perturbation theory assumption of weak coupling. In addition, a rate calculation assumes a particular form of kinetics, typically exponential, while quantum dynamics starts as a Gaussian (or cosine, e.g. Rabi oscillation) and may develop complex time-dependence, especially in the presence of multiple competing processes, such as charge- charge and charge-phonon scattering, charge and energy transfer, polaron and exciton formation, spin transitions, etc. that can occur on similar timescales. 6, 156-157 Further, Coulomb interactions and charge density distributions can exhibit strong dependence on atomic configuration, both in inorganic semiconductors 158 and especially in soft materials, such as metal-halide perovskites 18, 159 and organic matter. 160 Two approaches have been used to model the Auger-type quantum dynamics at the atomistic level of description and in the time-domain. One technique added intermediate high- energy states and decomposed a two-particle Auger transition into two single-particle transitions, as introduced by Kim et al. 156, 161 Zhu et al. used this method to demonstrate Auger-assisted electron transfer between semiconductor quantum dots (QDs) and molecular acceptors using non- adiabatic (NA) molecular dynamics (MD). 162 Wang et al. 13 introduced the global flux surface hopping (GFSH) method to allow NA transitions involving multiple particles simultaneously. Still, 47 both methods rely on single-particle coupling matrix elements and require many intermediate states in the calculation. The approaches have been used to model several systems, and proved capable of reproducing experimental timescales. 6, 104, 151, 163 Even though single-particle models have experienced huge successes, as exemplified by density functional theory (DFT) that maps a full many-body problem onto an effective single-particle description, 164 the physics of Auger processes resides in carrier-carrier interactions, and therefore, it is desirable to include such interactions explicitly into NAMD simulations. In this letter, we introduce the electron-electron interaction into NAMD, as implemented with ab initio real-time time-dependent DFT (TDDFT) in the PYXAID code. 15, 119 We compute the many-body Coulomb matrix elements and incorporate them into the TDDFT equations-of- motion and surface hopping (SH) transition probabilities. These matrix elements are readily available in linear response TDDFT codes. In contrast to linear response TDDFT, we avoid diagonalization of large matrices, keeping the computational cost low. By incorporating both NA and Coulomb coupling, the methods allow for simultaneous modeling of the quantum dynamic of electron-phonon and electron-electron scattering in a non-perturbative, configuration dependent way, in the time-domain and at the atomistic level of detail. We test the developed method on a CdSe quantum dot (QD), which has been studied for light harvesting and optoelectronics both experimentally and theoretically, and which constitutes a classic example of an Auger-type process in a nanomaterial. 142, 144, 150-151, 165-166 With the spatial confinement, QDs can exhibit large splittings of electronic state energies and are expected to have slow electronic transitions due to mismatch between electronic and vibrational energy quanta. However, most QDs exhibit rapid bulk-like relaxation. Our results show a sub-picosecond Auger process that outcompetes phonon- driven relaxation, in excellent agreement with experiments. Accurate and efficient, the developed 48 method can be applied to various low dimensional materials, in which charge-charge scattering provides an important pathway of energy flow. 4.2 Computational Details The evolution of the multi-electron wavefunction 𝚿 (𝒓 ,𝑹 ,𝑡 ) is determined by the time- dependent Schrödinger equation (TDSE): 𝑖ℏ 𝜕 𝚿 (𝒓 ,𝑹 ,t) 𝜕𝑡 =𝑯 (𝒓 ,𝑹 ,𝑡 )𝚿 (𝒓 ,𝑹 ,t) (4.1) The Hamiltonian H of the electron-nuclei system depends explicitly on electronic coordinates 𝒓 , and parametrically on the nuclear coordinates 𝑹 . 𝚿 (𝒓 ,𝑹 ,𝑡 ) can be represented with a linear combination of basis functions Φ 𝑖 (𝒓 |𝑹 ) : Ψ(𝐫 ,𝐑 ,t)=∑𝑐 𝑖 (𝑡 )Φ 𝑖 (𝐫 |𝐑 ) 𝑖 (4.2) In the time-dependent Kohn-Sham (TDKS) theory, 15, 118-119 the multi-electron wavefunctions Φ 𝑖 (𝒓 |𝑹 ) are Slater determinants of adiabatic Kohn-Sham 167 (KS) orbitals 𝜑 𝑘 that are computed with standard DFT packages, such as VASP 168 and Quantum ESPRESSO. 169 The TDSE is re- written in this representation as: 𝑖ℏ 𝜕 𝜕𝑡 𝑐 𝑖 (𝑡 )=∑(𝘀 𝑖 𝛿 𝑖𝑗 +𝑑 𝑖𝑗 )𝑐 𝑗 (𝑡 ) 𝑗 (4.3) where 𝛿 𝑖𝑗 is the Kronecker delta function, 𝘀 𝑖 is the eigen-energy for basis function Φ 𝑖 , and 𝑑 𝑖𝑗 is the non-adiabatic coupling (NAC) between basis states Φ 𝑖 and Φ 𝑗 : 𝑑 𝑖𝑗 =−𝑖ℏ<Φ 𝑖 | 𝜕 𝜕𝑡 |Φ 𝑗 >=−𝑖ℏ<Φ 𝑖 |∇ 𝐑 |Φ 𝑗 >𝐑 ̇ (4.4) Because 𝑑 𝑖𝑗 is a one-electron operator, it is non-zero only if Slater determinants Φ 𝑖 and Φ 𝑗 differ by one KS orbital. The NAC between multi-electron Slater determinants can be reduced to the NAC between single-particle KS orbitals. 15 49 In order to incorporate explicitly carrier-carrier interactions into the NAMD simulations, we introduce the corresponding many-body terms into the Hamiltonian 𝑯 (𝒓 ,𝑹 ,𝑡 ) . The electron- electron interaction between states Φ 𝑖 =Φ 𝑚 𝑛 =|𝜑 1 ...𝜑 𝑚 ...𝜑 𝑛 ⟩ and Φ 𝑗 =Φ 𝑝 𝑞 =|𝜑 1 ...𝜑 𝑝 ...𝜑 𝑞 ⟩ is evaluated as: 𝑉 𝑖𝑗 =⟨Φ 𝑚 𝑛 |𝑉̂ |Φ 𝑝 𝑞 ⟩=⟨𝑚𝑛 |𝑝𝑞 ⟩−⟨𝑚𝑛 |𝑞𝑝 ⟩ (4.5) Here, ⟨𝑚𝑛 |𝑝𝑞 ⟩ is the notation for the Coulomb integral: ⟨𝑚𝑛 |𝑝𝑞 ⟩=∫𝑑 𝒓 1 𝑑 𝒓 2 𝜑 𝑚 ∗ (𝒓 1 )𝜑 𝑛 ∗ (𝒓 2 )𝑟 12 −1 𝜑 𝑝 (𝒓 1 )𝜑 𝑞 (𝒓 2 ) (4.6) In particular, the matrix elements in Eq. (4.5) are non-zero for multi-particle states that differ in two electronic orbitals, and therefore, they describe electron-electron scattering. More complex types of Coulomb interactions involving three or more particles can be introduced using Bethe- Salpeter theory. 170 With the addition of the Coulomb matrix elements, the Hamiltonian 𝐻 𝑖𝑗 in Eq. (4.3) is updated to: 𝐻 𝑖𝑗 =𝘀 𝑖 𝛿 𝑖𝑗 +𝑑 𝑖𝑗 +𝑉 𝑖𝑗 (4.7) The NAC 𝑑 𝑖𝑗 , Eq. (4.4), describes electron-vibrational interactions, while the Coulomb matrix element 𝑉 𝑖𝑗 , Eq. (4.5), generates electron-electron scattering. Thus, both types of interactions are included into the TDSE. Hamiltonians of this type, Eq. (4.7), but without 𝑑 𝑖𝑗 , are commonly used in linear response TDDFT with hybrid functionals and many-body wavefunction models. 164, 170-174 In such cases, excited state energies are obtained as eigen-energies of the many-body Hamiltonian, and the wavefunctions are expressed as superpositions of many Slater determinants. Because of the large dimensionality of the many-body Hamiltonians, solution of the eigenvalue problem constitutes a 50 challenging problem, limited to small systems and few excited states. Our approach does not require finding eigenstates of the many-body Hamiltonian. Instead, the Hamiltonian, Eq. (4.7), is used to propagate the TDSE in the Slater determinant basis. An explicit consideration of a Slater determinant basis is not needed for solving real-time TDDFT equations. They can be solved numerically on a grid in the position or momentum space. 175-176 KS orbitals are used in real-time TDDFT as an auxiliary construct to deduce the electron density that is the key quantity of DFT. Adiabatic KS orbitals provide an alternative to the grid representation for solving the TDKS equations. The KS representation allows a great reduction of dimensionality, since a small set of relevant adiabatic KS orbitals can be selected by energy. 6, 104, 151, 163 The KS representation of the TDKS equations shifts the computational cost to solving the time-independent DFT problem. The main challenge resides in computing the NAC matrix elements, which can diverge in the case of trivial energy crossings. 177-178 A grid representation of the TDKS equations allows one to couple quantum mechanical electrons to classical nuclei in a mean-field manner, leading to the quantum-classical Ehrenfest approximation. 175, 179 The advantages and limitations of the Ehrenfest approximation are well documented. 12, 25, 116 In particular, the Ehrenfest approach works well for short-time dynamics. Simulating long-time dynamics of charge carriers in nanoscale systems requires inclusion of detailed balance between transitions upward and downward in energy. Detailed balance is needed for proper description of energy relaxation that leads to thermal equilibrium. Phonon-induced decoherence in the electronic subsystem should be included as well, since it can strongly influence transition times. Incorporation of the detailed balance and decoherence effects requires a choice of a basis, which is provided in our multi-particle simulations by the Slater determinants. While 51 these effects can be incorporated into the Ehrenfest method, 12, 116 we focus on the SH framework, which is commonly used for simulating NAMD The Hamiltonian, Eq. (4.7), can be used in any SH technique. We have implement the method within Tully’s fewest switches SH (FSSH) 10 that is the most popular SH approach, GFSH 13 that generalizes FSSH to allow super-exchange and direct many-particle transitions, and decoherence induced SH (DISH). 11 In particular, the transition probability for Tully’s FSSH is given by: 𝑃 𝑖 →𝑗 (𝑡 ,𝑑𝑡 )=∫ 2 ‖𝑐 𝑖 (𝑡 )‖ 2 𝑅𝑒 [( 𝑖(𝑑 𝑖𝑗 +𝑉 𝑖𝑗 ) ℏ )𝑐 𝑖 ∗ (𝑡 )𝑐 𝑗 (𝑡 )]𝑑𝑡 𝑡 +𝑑𝑡 𝑡 (4.8) where, 𝑑 𝑖𝑗 and 𝑉 𝑖𝑗 are the NAC and Coulomb matrix elements, Eqs. (4.4) and (4.5), and the coefficients 𝑐 𝑖 (𝑡 ) are solutions of the TDSE, Eq. (4.1), in the Slater determinant basis, Eq. (4.2). The simulations reported below are performed with the DISH approach, 11 which is rooted in the theory of quantum open systems, 115, 180 and naturally incorporates decoherence effects. Compared to FSSH and many other SH techniques that employ different ad hoc expressions for state-to-state transition rates, the transition probabilities in DISH are given by the quantum mechanical probabilities, i.e. squares of the 𝑐 𝑖 (𝑡 ) coefficients, directly. The approach is implemented within the PYXAID software package. 15, 119 Because KS orbitals 𝜑 𝑘 obtained from standard DFT packages are specified up to an arbitrary phase that can change between MD timesteps, we enforce phase-consistency. 181 In order to eliminate trivial crossings, 177 we keep track of orbital ordering, 178 although trivial crossings do not present a problem in the application considered here, because the KS orbitals are delocalized over the whole system. The simulations are performed with VASP using the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional and the projector augmented wave (PAW) method. 120-122, 168, 182- 183 The NAC are computed using the methodology developed recently for the PAW method 184 52 based on the Pawpyseed package. 185 The Coulomb matrix elements are also computed using the rigorous PAW methodology. 186 4.3 Case Study on CdSe Quantum Dot In order illustrate the advantages of incorporating the many-body Coulomb matrix elements into the NAMD simulations, we investigate electron cooling in the Cd33Se33 QD, Figure 4.1a. 10 Å of vacuum space is added along each direction to remove spurious interactions between QD images arising due to periodical boundary conditions. Only the Γ-point is used to sample the Brillouin zone, because the QD is a finite system. After heating and equilibrating the system at 300 K for 5 ps, a 2 ps MD trajectory is generated with a 1 fs timestep using the microcanonical ensemble. The NAC and Coulomb couplings are computed along the trajectory. The 2 ps trajectory is sufficient for simulating the Auger process, which occurs on a sub-picosecond timescale. The dynamics are fast, because the Coulomb interaction couples iso-energetic initial and final electronic states: the energy is transferred from electron to hole. For comparison, we performed simulations with just the NAC, turning off the Coulomb interaction. In such case, the electron energy can only be transferred to phonons, and the dynamics involving a transition across a large energy gap are slow. We also simulated intraband electron and hole relaxation, and electron-hole recombination. In the absence of a third charge carrier, the nonradiative recombination is slow and occurs by transfer of electronic energy to phonons. In order to simulate the slower processes, we replicated the 2 ps Hamiltonian multiple times under the classical path approximation. 15, 119 We illustrate the developed method with the CdSe QD, since it constitutes a classic example in which Auger-type energy exchange between electrons and holes breaks the bottleneck to the phonon-mediated relaxation. 142 The simulated Cd33Se33 QD has a diameter of around 1.4 nm, and its optimized structure is shown in Fig. 4.1(a). The electronic density of states (DOS), 53 computed for the optimized geometry, exhibits an isolated state at the conduction band (CB) edge, Fig. 4.1(b). The state, commonly denoted as 1Se, corresponds to the lowest unoccupied molecular orbital (LUMO) and couples weakly to other states. Fig. 4.1(c) demonstrates phonon-driven evolution of the energies of the highest occupied molecular orbital (HOMO), and the 1Se and 1Pe (LUMO+1) states along the MD trajectory. The 1Se state remains separated from the rest of the manifold at room temperature as well. The large gap (0.43 eV) between 1Se and 1Pe prohibits rapid non-radiative relaxation of hot electrons to the band edge (1Se), i.e. there exists a phonon- bottleneck for the hot electron relaxation. 139, 149 The bottleneck can only be observed experimentally, if the valence band (VB) hole is decoupled from the electron. 149 If both electron and hole are present in the QD and can exchange energy, the bottleneck is broken. 142, 151 The Auger channel is depicted schematically in Fig. 4.1(d). A hot electron deep inside the CB relax rapidly through the dense manifold of states to 1Pe by coupling to phonons. Then, the electron hops from 1Pe to 1Se, depositing its energy to the hole that is excited from the VB edge deeper into the band. After that, the hole relaxes rapidly to the VB edge through phonons. Finally, the electron and hole recombine across the bandgap, depositing the excess energy to vibrations. In the presence of more than a pair of charge carriers, e.g. a trion 187 or a double exciton, 155 the energy released during the recombination can be deposited into the extra carriers. These types of Auger processes are not considered in present model. 54 Figure 4.1. (a) Optimized structure of the Cd33Se33 QD (Cd: magenta, Se: green). (b) Electronic density of states for the QD. The DOS is higher in the valence than conduction band. The LUMO (1Se) is isolated from the LUMO+1 (1Pe) and the rest of the band. (c) Phonon-induced fluctuations of the orbital energies along the MD trajectory. (d) Schematic of the Auger-type process, in which the hot electron cools from 1Pe to 1Se while exciting the hole to deeper levels. The electron dynamics shown in Figure 4.2 demonstrate the sub-picosecond Auger process occurring in the system in the presence of the Coulomb interaction. Figure 4.2 shows both the multi-electron state populations and the corresponding single-electron state populations. The NAMD simulation is started with the excitation Φ 𝐻 𝐿 1 : HOMO → 1Pe. Fig. 4.2(a-b) presents the population evolution in the NAMD simulation with the Coulomb interaction. As shown in Fig. 4.2(a), state Φ 𝐻 𝐿 1 decays quickly, within ~0.2 ps, and at the same time, the group of states Φ 𝐷𝐻𝑠 𝐿 : (HOMO-1 and below → LUMO) gets populated. The evolution of the single particle populations shown in Fig. 4.2(b) further illustrates the Auger process. Within the first 0.2 ps of the simulation, the hot electron quickly jumps from 1Pe into the low energy level 1Se, and at the same time, the hole residing initially in the VB maximum, e.g. HOMO, gets excited into deep hole levels (DHs). Subsequently, the deep hole relaxes quickly, on a sub-picosecond timescale, back to the HOMO, as can be seen in the bell-shape curves in Fig. 4.2(a) for the population of Φ 𝐷𝐻𝑠 𝐿 and in Fig. 4.2(b) for DHs, and the inverse, V-shape curve for the HOMO. The Auger process followed by the quick hole relaxation bypasses the phono-bottleneck and accelerates the hot electron relaxation. 55 In comparison, the NAMD simulation done without the many-body electron interaction misses the Auger process, as evidenced by Fig. 4.2(c-d). In Fig. 4.2(c), the population of the states Φ 𝐻𝑠 𝐿 1 : (HOMO and below → LUMO + 1) decays slightly after 40 ps, and the hole excited states Φ 𝐻𝑠 𝐿 : (HOMO and below → LUMO) are not involved in the carrier dynamics at all. In Fig. 4.2(d), the population of the hot electron, initially residing in 1Pe, decays only slightly after 40 ps. On the other hand, the hole, which is initially in the HOMO, hops and equilibrates among the levels within kTB around the HOMO due to strong NAC between these states. Following the rapid initial equilibration of the hole near the VB edge, the DHs never get significantly populated, compare to Fig. 4.2(b). Without the many-body electron interaction, the electron decays from the 1Pe state to the 1Se state within hundreds of picoseconds, as observed in the experiments with the hole decoupled from the electron. 149 Therefore, the many-body electron interaction is essential for modeling of the Auger process. Figure 4.2. State populations illustrating the Auger process. (a), (c) Populations of the multi- electron states during NAMD with and without Coulomb scattering, respectively. Φ H L represents 56 the HOMO → LUMO excitation (blue). Φ Hs L1 represents all excitations from HOMO and below to LUMO+1, green). Φ DHs L represents all excitations from HOMO-1 and below to LUMO (red). (b), (d) Corresponding single particle populations during NAMD with and without Coulomb scattering, respectively. DHs (deep holes) represents all holes in HOMO-1 and lower states. Further insights into the mechanism of the Auger process and electron-hole relaxation are provided by Figure 4.3 and Table 4.1, which report the NAC values dij and the Coulomb matrix elements Vij. Figure 4.3 shows the NAC values between orbitals that may get populated in the NAMD simulation. Notably, the NACs between the HOMO, 1Se, and 1Pe are 1-2 orders of magnitude smaller than the NAC between states in the dense manifolds of the valence and conduction bands. The NAC for adjacent states is around 20 meV in the VB, and around 10 meV in the CB, rationalizing the fast hot-hole relaxation, and a slower hot electron relaxation, Figures S4.1 and S4.2. The NACs between the valence and conduction band edges are 100-fold smaller, less than 0.1 meV, as shown in Table 4.1 for the NAC between the ground state Φ 0 and excited states Φ H L and Φ H L1 . This is usually the case for wide gap systems and results in slow charge recombination. Here, the CdSe QD has a gap of 1.4 eV between HOMO and LUMO, and 0.43 eV between LUMO and LUMO+1. The NAC between 1Se and 1Pe, i.e. between Φ H L and Φ H L1 , is around 1 meV, much smaller than the value for the intra-band states. This weak phonon-electron scattering leads to the phonon-bottleneck that is expected to slow down the hot electron relaxation, and enhance its lifetime and the optoelectronic performance of QDs. However, the Coulomb matrix elements responsible for the electron-hole interactions are quite large, 20-40 meV between Φ H L1 and Φ Hn L (n=1, …, 14), Table 4.1. The Coulomb coupling provides an additional fast channel for the hot electron relaxation. The Coulomb matrix elements are on the same order or even larger than the NAC in the VB, rationalizing why the Auger process is faster the intraband hole relaxation, Figure 4.2(a),(b). 57 Figure 4.3. NAC between valence and conduction band orbitals. The orbital index is 0 for HOMO, -5 for HOMO-5, 1 for LUMO, 2 for LUMO+1, etc. The color and size of the squares show the magnitude of NAC between the corresponding orbitals. State I State J Coulomb Vij (meV) Φ H L1 Φ H1 L 40 ± 30 Φ H L1 Φ H14 L ~Φ H2 L (10~20) ± (10~20) State I State J NAC dij (meV) Φ 0 Φ H L 0.08 ± 0.06 Φ H L Φ H L1 1.5 ± 1.4 Φ 0 Φ H L1 0.06 ± 0.04 Table 4.1. Canonically averaged Coulomb matrix elements and absolute values of NAC matrix elements between selected states. Φ 0 is the ground state, Φ Hp Lq is the excitation from HOMO-p to LUMO+q. Figure 4.4 shows decay of the total electronic energy with and without the Coulomb scattering starting from the HOMO → LUMO+1 initial excitation. Figures S4.1 and S4.2 of Supporting Information show decays of energies hot electron and hole starting deeper in the bands. Table 4.2 summarizes timescales of the various processes. The electronic energy decay exhibits a 58 fast initial component with the Coulomb interaction turned on, Fig. 4.4(a). The electron exchanges energy with the hole, and the hole rapidly deposits the energy into phonons. Subsequently, a slow charge recombination follows. The Auger-assisted decay of the electron from 1Pe to 1Se occurs on a 0.15 ps timescale. The electron-hole energy exchange is fast because of the strong quantum confinement in the studied small QD. The Auger energy exchange is slower in larger QDs. 188 The energy transferred to the hole relaxes into phonons within 0.2-0.3 ps. Supported further by the charge-phonon relaxation data shown in Figures S4.1 and S4.2, this timescale is typical of intraband charge-phonon relaxation in nanoscale materials. 109, 189-190 The simulation without the Coulomb coupling does not exhibit the fast initial energy decay, Fig. 4.4(b). Instead, the 1Pe → 1Se transition of the electron requires 210 ps, representing the phonon bottleneck, characteristic of decoupled electron and hole. 139, 149 The electron-hole recombination occurs on a nanosecond timescale by coupling to phonons. An Auger-assisted electron-hole recombination requires another charge carrier that can accommodate the energy released during the recombination. This process in not considered in the present simulation. 59 Figure 4.4. Decay of the total electronic energy starting from the HOMO → LUMO+1 excitation (a) with and (b) without Coulomb scattering. The color strip shows the energy distribution during the NAMD simulation, and the dashed blue line shows the average energy. The ground state is used as the energy reference. With Coulomb Vij Without Coulomb Electron Relaxing to LUMO+1 - 0.3 ps Hole Relaxing to HOMO - 0.2~0.3 ps LUMO+1 → LUMO 0.15 ps 210 ps Electron-Hole Recombination 1000 ps 1000 ps Table 4.2. Timescales for hot electron relaxion to LUMO+1, hot hole relaxion to HOMO, LUMO+1 → LUMO (1Pe → 1Se) transition, and charge recombination. The hot electron/hole relaxations were modeled as single particle processes, not allowing energy exchange between particles. Simulations of the 1Pe → 1Se transition and electron-hole recombination included both electron and hole. Further insights into the Auger process are provided in Figure 4.5 that shows the evolution of the electron and hole energies, part (a), and the phonon influence spectrum for the intraband hole relaxation, part (b). In the first 0.1 ps, most of the energy lost by the hot electron is transferred to the hole, Fig. 4.5(a). However, a small fraction of the energy, several tens of meV, is already taken out of the electronic system, because the hot hole immediately starts transferring the energy to phonons. In general, the Auger and phonon driven processes proceed on similar timescales and compete with each other. 6, 156-157 The maximum energy reached by the hot hole is less than 1/3 of the energy released by the electron. The phonon influence spectrum shown in Fig. 4.5(b) indicates that the hole relaxes by coupling to optical phonons in the 20-40 meV frequency range (150-300 cm -1 ). The spectrum is computed by Fourier transforming the phonon-induced fluctuations of the energy gaps between HOMO and HOMO-n, where n=1…14, and averaging these fourteen Fourier transforms. 60 Figure 4.5. (a) Energy transfer from hot electron to hole during the Auger process. The red dashed curve shows the energy lost by the electron. The blue line shows the energy of the hole, which gains energy from the electron and then loses it to phonons. (b) Phonon influence spectrum for hot hole relaxation, averaged over transitions between Φ H L1 and Φ Hn L : (HOMO-n → HOMO, n=14, 13, …, 1). 4.4 Conclusions In summary, we have developed and demonstrated an efficient technique for modeling Auger-type processes with NAMD. The quantum dynamics calculation incorporates simultaneously electron-phonon and electron-electron scattering processes in a non-perturbative manner with explicit dependence on atomic configuration. In addition to the NAC that describes electron-vibrational interactions, we have incorporated the Coulomb coupling that captures many- body electron interactions. In particular, the Coulomb matrix elements can couple states that differ by two KS orbitals, and therefore, allow transitions in which two particles exchange energy. The methodology has been implemented within several SH schemes. The additional expense of the calculation arises from the need to compute configuration dependent Coulomb matrix elements, which are computed routinely in linear response TDDFT. In contrast to linear response TDDFT, the current method does not require diagonalization of large matrices, and the additional expense is minimal. 61 The developed technique has been demonstrated by application to the Auger-assisted relaxation of hot electrons in a CdSe QD. This problem constitutes a classic example of an Auger process in a nanoscale material. The Auger-type energy exchange between electron and hole breaks the phonon-bottleneck to the hot electron relaxation. The bottleneck can be recovered only when the hole is decoupled from the electron. Our simulations demonstrate all processes observed experimentally, and reproduce the experimental timescales. Importantly, Auger- and phonon- assisted processes can occur in parallel on similar timescales and compete with each other. Therefore, it is necessary to model them simultaneously, as achieved by the developed technique. The method makes no perturbative assumptions, as is often the case in rate calculations. This is important, because Coulomb interactions can be strong in confined systems. Further, the quantum dynamics simulations enabled by the developed approach differ from kinetics calculations that have to assume a particular kinetic mechanism, e.g. 1 st or 2 nd order, Gaussian or exponential. The NAC and Coulomb terms depend explicitly on nuclear configuration, which is particularly important for studying softer materials, such as metal halide perovskites and organic matter, as well as in organic materials with surfaces, edges and defects, since such systems can undergo large-scale anharmonic motions as a result of both a finite temperature and an electronic excitation. 191-192 The direct modeling of Auger processes with NAMD captures realistic aspects of materials’ atomic structure and provides an efficient and accurate approach for studying carrier dynamics in low dimensional materials. 62 Chapter 5: Unsupervised Machine Learning of Nonadiabatic Molecular Dynamics in MAPbI3 Reprinted with permission from ACS Ener. Lett. 2020, 5, 6, 1930-1938, Copyright 2020, American Chemical Society. Authors: Guoqing Zhou, Weibin Chu, Oleg V. Prezhdo Contributions: O.V.P. designed the project, G.Z. implemented the model and analyzed the results, W.C. performed the simulations. 5.1 Introduction Hybrid organic-inorganic perovskite solar cells (PSC) have attracted significant attention over the last decade because of their extraordinary properties, such as high power conversion efficiency, and ease and low cost of production. 193-198 Halide perovskites is a promising material for the next generation photovoltaic devices due to long electron and hole lifetimes, large charge diffusion lengths, low defect densities, and high absorption coefficients. 198-200 Applications of perovskites extend beyond solar energy into light-emitting diodes. 201-206 The efficiency of PSC has reached over 23%. 193, 207 Tremendous efforts continue to be dedicated to enhancing the efficiency, improving stability and reducing environmental impact, by employing different synthesis methods, 208-210 adapting new device architectures, 211-212 and trying various stoichiometries. 202, 213- 214 Multiple experiments have been dedicated to understanding the unique combination of properties that make PSC so efficient. 198, 215-218 The unusual charge carrier dynamics and its response to nuclear motions are among the main factors determining the PSC performance. Multiple theoretical efforts have focused on elucidating the geometric and electronic structure of perovskites, focusing on different stoichiometries, dimensionalities (3D vs. 2D vs. quantum dots), 63 and defects (vacancies, interstitials, dopants, grain boundaries, interfaces). 125, 199, 213, 219-221 Theoretical studies of the quantum dynamics of coupled electronic and nuclear degrees of freedom are still rare, since they require advanced computational techniques. 15, 26, 118-119 Analysis and interpretation of simulation results constitute an additional challenge, in particular, due to the complex landscape of perovskite compositions, structures, length- and time-scales involved. Machine learning (ML) has emerged as a useful tool for analyzing both large and limited amounts of data available from either experiment or modeling, with applications in an extremely diverge range of fields, from high energy physics, to chemistry and materials science, to social sciences, politics, etc. 222-225 Often, ML allows scientists to uncover trends and correlations that are hard to anticipate a priori. Quantum dynamics of charge carrier trapping and recombination in perovskites is an excellent candidate for ML analysis, because of the multitude of factors influencing carrier lifetimes, including perovskite composition, thermodynamic phase, defects, geometric and electronic structure, temperature, atomic motions, electron-phonon coupling, quantum coherence, and so on. Generally, ML techniques can be categorized into supervised and unsupervised learning. 225 Computational studies of nanoscale materials, and perovskites in particular, can benefit from both approaches. For example, supervised learning can be used to obtain efficient neural-network force-fields trained to ab initio density functional theory (DFT), 16, 226 while unsupervised learning can uncover complex reaction paths. 227 ML techniques help resolve the complex dependence problem and attract strong attention from material scientists, since they provide analysis and modeling tools that outperform the traditional methods in terms of accuracy and efficiency. Nonadiabatic molecular dynamics (MD) simulation combined with time-dependent DFT (TDDFT) has recently established itself as a tool of choice for modeling the complex highly non- 64 equilibrium evolution of coupled electronic and nuclear degrees of freedom in molecular, nanoscale and condensed matter systems. 15, 119, 135, 228-230 Nonadiabatic MD treats all particles and interactions explicitly and at the atomistic level, and avoids approximations used in phenomenological theories, such as harmonic phonons, weak system-bath interaction, few level models, as well as adjustable parameters. By generating a time-domain description of excited state quantum dynamics Nonadiabatic MD - TDDFT mimics directly multiple time-resolved pump- probe spectroscopies, while at the same time taking explicit account of the realistic aspects of material structure, including chemical composition, defects, dopants, surfaces, interfaces, grain boundaries, etc. The ab initio time-domain simulations are easily transferrable to new types of materials. Since such simulations produces large amounts of data, it is natural to consider ML as a data analysis tool. In this letter, we apply unsupervised ML to ab initio Nonadiabatic MD and demonstrate that nonradiative charge carrier losses in MAPbI3 are governed primarily by material’s structural properties rather than particular phonon motions. This fact is unexpected because the nonadiabatic electron-phonon coupling depends explicitly on nuclear velocity. Also unexpectedly, we find that MA motions have a strong influence on the nonradiative relaxation, even though MA does not contribute to the electron and hole wavefunctions. The main influence arises from geometry of the inorganic Pb-I sublattice, with the I-I-I angle providing a stronger measure of electron-phonon coupling and bandgap than the Pb-I-Pb tilt angle. These findings rationalize the unusual temperature dependence of charge carrier lifetimes in halide perovskites 231 and support the mechanism of polaron formation involving Pb-I lattice distortion and MA rotation. Our study highlights the ability of ML to uncover important correlations that are not anticipated a priori within complex datasets, and shows that it can provide important quantitative measures with 65 limited amounts of data. The reported findings suggest that PSC performance can be improved by engineering of perovskite geometric structure using site substitution, doping, mixed stoichiometries and related synthetic techniques. 5.2 Computational Details 5.2.1 Non-Adiabatic Molecular Dynamics Figure 5.1. Geometric and electronic structure of MAPbI3. (a) Tetragonal phase of MAPbI3 (tetragonal primitive unit cell) (b) Projected density of state showing inorganic, Pb and I, and organic CH3NH3 (MA) contributions. (c) HOMO and (d) LUMO charge densities, localized primarily on I and Pb, respectively. 66 We study pristine tetragonal MAPbI3 represented by one 48-atom 111 supercell as shown in Figure 5.1a and a larger 384-atom 222 supercell to check the convergence of the results. The 1x1x1tetragonal phase is represented using a √2√22 unit cell of the cubic phase. DFT and ground state MD simulations are performed with the Vienna Ab initio Simulation Package (VASP), 120, 168, 182-183 employing the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional 121 and projector augmented wave (PAW) method. 122 We use 4×4×2 Monkhorst−Pack k-point mesh to sample the Brillouin zone for structure optimization, static calculation and MD simulation. For the 2x2x2 supercell we only use gamma due to the computational cost. The optimized structure has bandgap 1.63 eV and lattice constants a = 8.8 Å, c = 12.685 Å, in good agreement with the previous publications and experimental values. 232-233 Since MAPbI3 is a direct bandgap semiconductor and we mainly focus on the nonadiabatic recombination, only the wavefunction from -point is used to compute the Nonadiabatic Coupling (NAC) for the machine learning analysis. The crystal is first optimized and then thermalized at 300 K by velocity rescaling. Then, a 9 ps ab initio MD trajectory is generated with a 1 fs timestep, and the NAC between HOMO and LUMO -point is computed using the PYthon eXtension for Ab Initio Dynamics (PYXAID) code. 15, 119 The PYXAID methodology has been used previously to model successfully excited dynamics in a broad range of materials 104, 109, 124, 126, 132-133, 234-236 , the results of the quantum dynamics calculations performed with PYXAID on halide perovskites also show good agreement with the experimental data. 237-238 Generally, the timescales and mechanisms of charge carrier dynamics depend on the energy gap, NAC and coherence time between pairs of states. The latter can be computed from autocorrelation of the energy gap. 123 Therefore, we focus on the gap and the NAC, and apply ML to analyze how they depend on atomistic features and motions of MAPbI3. 67 5.2.2 Mutual Information A total of 44 features are considered for both the 111 supercell and the 222 supercell, including 22 static ones and 22 dynamic ones, which are the time derivative of the static ones, as listed in the Table S5.1 and Fig. S5.1 in the supporting information. The 22 static features, like bonds, bond angles, dihedral angles, relative orientations and distances, are represented by 207 specific internal coordinates in this 48-atom system. For example, as for the feature C-N bond distance, there are 4 C-N bonds from the 4 MA molecules. As a result, 414 variables are extracted from the MD trajectory. The top 10 features are selected and discussed in detail. The ML analysis demonstrates that the importance of the structural and dynamic parameters drops rapidly within the set of 44 features, with the 10 th most important feature already being 1.5-1.6 times less important than the top feature. Figure 5.2 shows distributions of the top 3 features (bond angles within the inorganic sublattice of MAPbI3) and the target variables (NAC and bandgap). The 414 variables and the target variables obtained from the MD trajectory are normalized individually to remove bias and artificial variance introduced by the choice of units. After that, the pairwise mutual information (MI), I(X,Y)=∬d𝑥 d𝑦 𝑝 (𝑥 ,𝑦 )log ( 𝑝 (𝑥 ,𝑦 ) 𝑝 (𝑥 )𝑝 (𝑦 ) ) , (5.1) between each feature and each target variable is estimated using K-nearest neighbor (KNN) algorithm. 239-240 In addition, the entropy is estimated as H(𝑋 )= −∫d𝑥𝑝 (𝑥 )log (𝑝 (𝑥 )) , (5.2) and the MI is normalized by dividing by the entropy of the corresponding target variable. In the above equations, p(x,y) is the joint density distribution of continuous random variables X, Y, and p(x), p(y) are their densities. Base 2 is used for the logarithm calculation. The results of the KNN algorithm calculations reported below use k = 3. We have tested different k values and observed 68 similar results. A smaller k leads to a shaper fluctuant probability distribution, but introduce noise/variance to the estimated MI. A larger k gives smoother densities; however, it is harder to differentiate contributions of different variables. With the computed MI values for these 414 variables, averaging is done across the variables for each type of features, to get the final MI values for each feature. 5.3 Results and Discussion Figure 5.1 characterizes the geometric and electronic structure of MAPbI3 in its optimized geometry. The MAPbI3 bandgap falls within the solar spectrum, which is a necessary condition for solar energy applications. Spatial separation of electrons and holes is one of the reasons for the fantastic performance of PSC. The valence band edge arises mostly from atomic orbitals of iodine atoms, while the conduction band edge in on lead atoms. The sparse structure of the inorganic Pb- I lattice, containing sufficiently large cavities to accommodate MA and similar molecules, ensures that the overlap of the atomic orbitals supporting electrons and holes is small, compared with the more conventional two-component semiconductors, such as TiO2 and CdSe, and pure elements, such as Si. The small overlap ensures small NAC and short quantum coherence time, favoring long-lived charge carriers. 213, 220, 231 Figure 5.2 presents canonical distributions of the NAC and the bandgap, as well as the three most important MAPbI3 features that share the highest MI with the NAC and the bandgap. These distributions are broad and roughly Gaussian, implying that thermal fluctuations of the lattice have a strong impact on the electronic system, motivating us to search for correlations of lattice coordinates and motions with the electronic properties. Considering each panel in Figure 5.2, we note that the NAC is indeed very small, usually less than 2 meV. However, the NAC has fairly long tails that may have a significant influence on charge carrier recombination, since the 69 recombination rate depends on NAC squared, for instance, according to Fermi’s golden rule. 241 The bandgap fluctuates by about 0.5 eV at 300 K. The fluctuation is large because the conduction and valence band edges originate from spatially separated and weakly interacting atomic orbitals, Figure 5.1. The large bandgap fluctuation gives rise to short quantum coherence for the electron- hole recombination. 123 A short coherence slows down quantum dynamics and extends carrier lifetimes. 139 The I-I-I and I-Pb-I angles form Gaussian distributions around 90, Figure 5.2c-d, with the former exhibiting a larger fluctuation. The Pb-I-Pb angle is notably less than 180. This octahedral tilt angle plays an important role in polaron formation. 242-243 Figure 5.2. Probability distributions of (a) NAC between HOMO and LUMO, meV, (b) bandgap, eV, and (c) I-I-I, (d) I-Pb-I, (e) Pb-I-Pb angles, degrees. These angles have the highest mutual information with NAC and bandgap among all considered geometric features and atomic motions. One can attempt to establish a relationship of the NAC and the bandgap with geometric features and motions of MAPbI3 using linear models. However, even though each individual 70 feature has a trivial Gaussian-like distribution, their joint distributions are complicated, and there is no obvious linear dependence between them. Figure 5.3 presents the time evolution and joint distributions for the three coordinates with the NAC. The NAC fluctuates faster than the coordinates, which vary more slowly and smoothly. The largest absolute value of the correlation coefficient between the geometric features and the NAC is less than 0.32, with the majority of the correlation coefficients around 0.1 to 0.2. The fairly computationally expensive 9 ps ab initio MD trajectory provides 9000 data points, which is sufficient to evaluate the probability distributions for individual properties, Figure 5.2, and to converge Nonadiabatic MD calculations. However, the joint distributions are distorted, Figure 5.3d-f, making it harder to determine the contributions from the atomic features to the target variables. Figure 5.3. Evolutions of (a) I-I-I, (b) I-Pb-I(90), (c) Pb-I-Pb angles superimposed with evolution of NAC between HOMO and LUMO. (d), (e), (f) Joint distributions of these three angles with NAC. Unsupervised learning techniques can extract hidden complex information from limited data, providing an important advantage over the traditional analyses. Tables 5.1 and 5.2 show the 71 most important internal coordinates and motions that have the highest MI with the NAC and the bandgap, respectively. The entropies of the NAC and the bandgap are 1.84 ± 0.04 and 1.98 ± 0.04, respectively, which are close to the value of the standard normal distribution log 2 (2𝜋𝑒 )/2 = 2.05. The MI is used to quantify the shared information between two random variables. The I-I-I angle shares the highest information with the NAC and the bandgap, MI = 0.87 ± 0.05 and 2.31 ± 0.09, respectively. The MI values are much larger than the MI of the bivariate normal distribution, MI= −log 2 (1−𝜌 2 )/2 < 0.07, calculated using the correlation coefficient |𝜌 | < 0.32 in our data set, Figure 5.3d. This result points to the complex and significant interdependence of these features, and demonstrates invalidity of linear correlation. Table 5.1. Mutual information of internal coordinates and motions with NAC I-I-I angle I-Pb-I (90) angle Pb-I-Pb angle I-Pb-I (180) angle Pb-I-Pb motion rank 1 2 3 4 5 MI 0.87±0.05 0.73±0.06 0.71±0.07 0.63±0.04 0.62±0.06 MI/H 0.48±0.03 0.39±0.04 0.39±0.04 0.34±0.03 0.34±0.04 MA-axis c angle I-Pb-I (90) motion MA-MA angle MA-MA distance MA-axis a angle rank 6 7 8 9 10 MI 0.61±0.08 0.60±0.05 0.60±0.03 0.59±0.02 0.59±0.04 MI/H 0.33±0.04 0.33±0.03 0.33±0.02 0.32±0.01 0.32±0.03 The top 10 features exhibiting the highest MI with the NAC, Table 5.1, include, in decreasing importance: I-I-I bond angle, I-Pb-I bond angle (~90), Pb-I-Pb bond angle, I-Pb-I bond angle (~180), angular velocity of the Pb-I-Pb bond angle, MA orientation relative to the lattice vector c, Figure 5.1, angular velocity of the I-Pb-I bond angle (~90), relative angle between two neighboring MA molecules, center-of-mass distance between two neighboring MA molecules, and orientation of MA relative to the lattice vector a. These 10 most important features are selected by the unsupervised learning from the set of 44 features, including internal coordinates and motions 72 of bond lengths, bond angles and dihedral angles in the inorganic Pb-I lattice and MA molecules, as well as orientations of MA and distances between them. The bond angles within the inorganic lattice have the highest MI among all features. Note that while the Pb-I-Pb tilt angle, frequently used in perovskite structure analysis, 242-243 belongs to the top 10 features, the I-I-I angle provides a better measure of the impact of lattice distortion on charge carrier dynamics. The tilt and I-I-I angles together describe the relative orientation of the PbI6 octahedra blocks, and the I-Pb-I angles (90 and 180) describe the PbI6 octahedra shape. This is consistent with the data in Figure 5.1. The HOMO and LUMO of the system are localized on the inorganic lattice. Therefore, the NAC between these orbitals is mostly determined by the lattice angles. The above results suggest that modifying the structure of the PbI6 octahedra block, for instance, by strain, temperature, chemical composition or doping, should influence the carrier dynamics. Perhaps most surprisingly, the ML analysis indicates that MAPbI3 lattice motions are much less important than the MAPbI3 structure. Only two motions are present among the top 10 features in Table 5.1. This result is unexpected, because the NAC depends explicitly on nuclear velocity 𝑅 ̇ : NAC=−𝑖ℏ<𝜙 𝑖 | 𝑑 𝑑𝑡 |𝜙 𝑗 >=−𝑖ℏ<𝜙 𝑖 |∇ 𝑅 |𝜙 𝑗 >∙𝑅 ̇ (5.3) This expression separates the influence of the velocity and overlap contributions. The ML data demonstrate that the overlap-type matrix element <𝜙 𝑖 |∇ 𝑅 |𝜙 𝑗 > is more important for the charge carrier dynamics in MAPbI3 than the nuclear velocity. This fact explains the unusual temperature dependence of the charge carrier lifetimes. 231, 244 The velocity dependence of the NAC suggests that charge carrier lifetimes should decrease with increasing temperature, since higher temperature implies higher kinetic energy and faster motions and larger NAC, leading to faster charger carrier relexation. Instead, the lifetime increases with temperature, 244 because thermal fluctuations induce lattice distortions and decrease <𝜙 𝑖 |∇ 𝑅 |𝜙 𝑗 >. 231 73 Also surprisingly, 4 out of the 10 top features involve MA molecules, Table 5.1, even though they do not contribute to the valence and conduction band edges, Figure 5.1b. These molecules influence the charge carrier dynamics in two ways. Carrying a total positive charge and having a non-uniform charge distribution, they interact with electrons and holes electrostatically. Rotations of MA molecules contribute strongly to polaron formation, 242-243, 245-247 screen electron- hole interactions, 248-249 and influence ion diffusion responsible for the current-voltage hesteresis. 250 The screening facilitates electron-hole separation, reduces charge carrier recombination and increases carrier lifetime. In addition to the long-range electrostatic effects, MA molecules couple to the inorganic lattice by short-range repulsive steric interactions, and affect electron-hole dynamics indirectly by distorting the inorganic lattice. The internal coordinates of the MA molecules contribute little to the carrier dynamics, and their MI with the NAC and the bandgap is small. Internal stretching and bending MA motions have little effect on both electrostatic and steric interactions. Table 5.2 presents the MI between the internal coordinates and the bandgap. The motions of the internal coordinates are not included into consideration, since they do not affect the band gap, in particular, since the DFT energy terms do not depend on nuclear velocities. Table 5.2 leads to similar conclusions as Table 5.1. The I-I-I angle shares the largest MI with the bandgap, just as with the NAC. The orientation and shape of the PbI6 block, as well as the orientation of MA molecules are the most important factors affecting the electronic properties of the system. The last entry in Table 5.2 is the Pb-I bond length, which does not enter Table 5.1. The MI values drop significantly after these 9 most important features, and we do not present them here. 74 Table 5.2: Mutual information of internal coordinates with bandgap I-I-I angle I-Pb-I (90) angle Pb-I-Pb angle MA-axis c angle I-Pb-I (180) angle rank 1 2 3 4 5 MI 2.31 ± 0.09 1.99 ± 0.11 1.96 ± 0.14 1.72 ± 0.16 1.69 ± 0.07 MI/H 1.17 ± 0.05 1.00 ± 0.06 0.99 ± 0.07 0.87 ± 0.08 0.85 ± 0.04 MA-MA angle MA-MA distance MA-axis a angle Pb-I distance rank 6 7 8 9 MI 1.67 ± 0.13 1.62 ± 0.06 1.56 ± 0.12 1.47 ± 0.09 MI/H 0.84 ± 0.07 0.82 ± 0.03 0.79 ± 0.06 0.74 ± 0.05 The top 10 features are almost the same from the results of the large 384-atom 222 sueprcell system and the small 48-atom 111 system. Most importantly, the first 3 features are the same for the two systems and for the two target variables, the NAC and bandgap. The top 7 out of 10 features sorted based on MI with NAC are all internal coordinates rather than motions. The same resutls confirm that I-I-I angle is a better structure parameter governing the nonadiabatic charge recombination in MAPbI3 and that the structure deformation is more important than the lattice variation regarding the charge recombination. However, one key difference is Pb-Pb, which is the distance between Pb from the first and second nearest neighbors. Pb-Pb is not in the top 10 list for the small system, but Pb-Pb is 4 th place for the correlation with NAC and bandgap, as show in the Table S5.3 and S5.5. This difference is due to the strong correlation between the nearby Pb atoms in the 111 system, which reduced the deformation of these octahedra blocks, specifically the distance between these blocks represented as Pb-Pb distance, and this results in less contribution of Pb-Pb to the NAC and bandgap in this system. Moreover, it is shown by Lahnsteiner et al that artificial correlation is induced to the reduced supercell size system in MAPbI3, 251 which will affect the MA motion, and also the Pb-Pb deformation as well. This shows 75 that the deformation of octahedra blocks, as represented by I-I-I, I-Pb-I, Pb-I-Pb angles, and Pb- Pb distance, play a key role in the nonadiabatic carrier relaxation in MAPbI3. Another traditional empirical way to evaluate the influence of various coordinates and motions on the NAC and the bandgap is go to the frequency domain. Figure 5.4 shows Fourier transforms (FTs) of the phonon-induced fluctuations of the NAC and the bandgap along the ab initio MD trajectory, superimposed with the FTs of the three main angles. The spectra overlap significantly, indicating further that the bandgap and the NAC are correlated to these three features. In all cases, the strongest FT signals arise in the low frequency region (< 100 cm -1 ). FT of the I-I- I angle fluctuation matches best the higher frequency part of the NAC and bandgap FTs, while the I-Pb-I and Pb-I-Pb angles provide a better match at the lower frequencies. The NAC exhibits stronger signals at the higher frequencies compared to the bandgap, because the NAC is a higher- order property than energy and is sensitive to a broader range of motions. Figure 5.4. NAC and bandgap spectra compared with spectra of the key internal coordinates: I-I- I, I-Pb-I (90) and Pb-I-Pb angles. 76 5.4 Conclusions In conclusion, we have reported a ML analysis of the NAC and the energy gap that determine nonradiative charge carrier lifetimes in MAPbI3. The study shows that MAPbI3 structure, rather than its motions, governs most strongly the electron-phonon interactions responsible for charge carrier losses. This result is surprising, because the NAC depends explicitly on nuclear velocity, and therefore, faster and larger amplitude motions should enhance the nonradiative relaxation. This result rationalizes, in particular, the unusual temperature dependence of charge carrier lifetimes in halide perovskites. Charge losses accelerate at higher temperatures in most materials, however, they slowdown in MAPbI3 because of increased structural deformations. Only 2 motions are present among the top 10 MAPbI3 features governing nonradiative charge recombination. The most important features stem from the inorganic Pb-I lattice as expected, since electrons and holes localize primarily on Pb and I atoms. Unexpectedly, 4 of the 10 features involve rotations and displacements MA species, which carry no contribution to the conduction and valence band edges of MAPbI3. This finding can be explained by direct long-range electrostatic interactions between charge carries and the MA cation, and indirectly by short-range steric interactions between MA and the Pb-I lattice. Internal motions of MA have little influence on the charges. ML provides a better tool than simple regression models for analyzing problems with complex dependencies. ML gives quantitative results that are hard to generate with the traditional methods such comparing spectral densities. ML performs well with limited amounts of data, when spectral densities and correlations show significant noise. This is particularly valuable with costly ab initio time-domain simulations. The noise is partially eliminated during the calculation of MI 77 by the summation or integration. The results of the unsupervised learning analysis with MI demonstrates the capability of extracting critical factors affecting carrier losses PSC. The fact that structural deformations of the inorganic lattice are most important among all dynamic and static features in PSC suggests that the PSC performance can be engineered by doping or site substitutions affecting the structure. The ML analysis shows that the I-I-I angle provides a better measure than the Pb-I-Pb tilt angle when studying the effects from these modifications on charge carrier lifetimes. 78 Chapter 6: GPU-Accelerated Semi-Empirical Born Oppenheimer Molecular Dynamics using PyTorch Reprinted with permission from J. Chem. Theory Comput. 2020, 16, 8, 4951-4962, Copyright 2020, American Chemical Society. Authors: Guoqing Zhou, Ben Nebgen, Nicholas Lubbers, Walter Malone, Anders MN Niklasson, Sergei Tretiak Contributions: B.N. and S.T. designed the project, G.Z. implemented the model and performed all the testing and analyzed results. A.M.N. provides the demo SP2 and XLBOMD codes. All contributed to the discussion and preparation of the manuscript. 6.1 Introduction Semi-Empirical Quantum Mechanics (SEQM) models are widely used to efficiently study ground and excited state electronic properties of chemical systems and materials. 252-253 While high level methods like configuration interaction (CI), coupled cluster (CC), and density functional theory (DFT) can have a higher accuracy than SEQM, the former techniques scale cubically to exponential with the system size and have a large computational overhead, typically limiting them to systems with 10-100 atoms. 30, 167 SEQM models utilize very minimal basis sets and parameterize various integrals in the electronic Hamiltonian with experimental or ab initio references, 254-258 resulting in methods that are much faster than traditional wavefunction or DFT frameworks while maintaining a reasonable level of accuracy. The increased speed of SEQM techniques allow them to be applied to systems traditionally outside the reach of quantum mechanical methods such as proteins, 253 nanotubes, 259 and many noncovalent complexes, 260-261 . However, despite the many decades of research into SEQM methods, 252, 254-255, 257-258, 262-263 the limitations of the functional 79 form for the underlining reduced electronic Hamiltonian models leave much room for improvement in their accuracy. 257, 263 Over last several decades, several distinct groups of semiempirical methods have been developed, including MNDO (Modified Neglect of Diatomic Overlap), 254 AM1 (Austin Model 1), 255 MNDO/d, 262 PM6 (Parametric Model 6) 257 and PM7. 263 These approaches were developed in series to improve the description of core-core interactions and to produce more accurate parameter schemes. To accurately describe non-covalent bond interaction, Hobza and co-workers extended PM6 with models like PM6-DH, 264 PM-DH2, 265 PM-DH+. 266 Moreover, a series of recently developed Orthogonalization-corrected methods (OMx), reach higher accuracy and have better descriptions of non-covalent bonds by taking into account repulsive orthogonalization, attractive penetration effects, repulsive core-valence and dispersion interactions. 258 While these developments have resulted in significant improvements in effective Hamiltonian models, they still lack the accuracy and generality compared to conventional high level quantum methods. 267- 268 This poses a need for further improvements of SEQM using modern data science approaches, as was demonstrated for the parameterization of reduced Density Functional Tight Binding (DFTB) models. 269 In the last couple of decades, graphics processing units (GPUs) have demonstrated their ability to accelerate numerically intense computations. The matrix operations that comprise many parts of SEQM methods can strongly benefit from such GPU offloading. 252 For example, the most time consuming component of SEQM for mean-field ground state calculations, the Self-Consistent Field (SCF) procedure, can be speed up significantly through the use of novel numerical methods, such as the single-particle density matrix expansion algorithm with second order spectral projection polynomials (SP2). 270 GPU implementations of these approaches, which used to require 80 highly specialized knowledge of the layout and structure of GPUs, has become much simpler through the use of high-level packages such as PyTorch, 271 a package originally designed for machine learning applications. Further, the incorporation of reverse-mode Automatic Differentiation, 271 which can compute the gradient of a quantity with respect to all parameters at the same computational cost as computing the quantity itself, has facilitated the solution of complex optimization problems. All of these developments have contributed to the rising popularity of machine learning in many areas including but not limited to physics, chemistry and material science. Applications in these fields range from predicting various properties including atomic partial charge, 29 molecular diploes, 17 atomization energy, 28 generating fast and highly accurate potential energy surfaces and forces, 16 modeling chemical reactions, 272 and discovering new materials. 273 Unfortunately, many of these applications suffer from the “black box” problem, where a deep neural network (DNN) makes an accurate prediction, but provides essentially no explanation as to why that prediction was made. A possible solution to this transparency problem is to build quantum mechanical models on top of a DNN for modeling molecular systems. Yaron et al. 269 demonstrated a DNN linked with a self-consistent DFTB layer and successfully reduced the error in DFTB by 40% to 60% when predicting total molecular energy and electronic dipole moment. Zubatyuk and coworkers linked a simple extended Hückel model with HIPNN 274 to dynamically generate parameters for the semi-empirical model which facilitated electronic structure predictions comparable to DFT calculations. 275 Future work similar to these developments would tremendously benefit from a general software platform combining an efficient implementation of model Hamiltonians in a programming platform which also natively supports machine learning techniques; SEQM aims to fill this role. 81 In this article, we present an implementation of various Semi-Empirical Quantum Mechanical methods utilizing PyTorch (PYSEQM) in an open source software package. 276 PYSEQM supports a variety of SEQM Hamiltonians including MNDO, AM1, PM3, PM6, and extended Hückel theory. 271 Additionally, it contains a simple Molecular Dynamics (MD) driver to facilitate the simulation of dynamic properties and non-equilibrium configurations for molecular systems. Furthermore, PYSEQM is capable of processing multiple molecules simultaneously on a GPU (i.e., run multiple MD trajectories in parallel) to get energy and gradient information, which may be of use in future machine learning (ML) applications or in various theoretical methodologies requiring ensemble averaging of a swarm of trajectories. 277 Since PyTorch supports a wide range of execution options on GPUs, PYSEQM can easily be deployed on a wide range of GPU devices to achieve high performance. Finally, we have implemented Extended Lagrangian Born Oppenheimer Molecular Dynamics (XL-BOMD) 278-282 for accelerated MD simulations that avoids the costly self-consistent field (SCF) algorithm as well as the recursive Fermi-operator density matrix expansion algorithm (SP2) for a rapid SCF convergence utilizing GPU devices. 270, 283-285 The largest limitations of PYSEQM are that it currently only works for closed shelled systems and that it only supports s and p-orbitals, though an extension to d-orbitals is currently under development. Additionally, the SP2 scheme currently implemented in PYSEQM will fail with bond breaking; adding newer versions of the scheme would address this shortcoming. 286 We detail the implementations of these semi empirical models as well as XL-BOMD and SP2 in Section 2. In Section 3, we show the performance of PYSEQM on two groups of systems: the nanostar dendrimer 287 and various lengths of polyethylene chains. The testing on the nanostar demonstrates the performance of PYSEQM on large systems in addition to the enforcement of energy conservation with the XL-BOMD scheme. The second test set with polyethylenes shows 82 the scaling of this code with system size and demonstrates situations where PYSEQM outperforms conventional semi-empirical codes. Finally, we summarize the results and conclude in Section 4. 6.2 Methods In this section, we detail the features of PYSEQM package. First, we examine the form of the three unique semi-empirical methods implemented in PYSEQM: MNDO, AM1 and PM3 for elements from H to Cl covering organic molecular systems. While not explored in this paper, PM6 with d-orbitals is also fully implemented. These effective Hamiltonians are constructed with a set of empirical parameters fit to experimental data to compensate for the neglect of electron- correlation inherent in all single reference (mean-field) methods. Second, we detail the SP2 method for rapid density matrix determination directly from a Hamiltonian matrix, bypassing matrix diagonalization. SP2 expands the density matrix in terms of the Hamiltonian Operator with an iterative scheme, which is much more efficient on modern GPU architectures than traditional diagonalization and density matrix construction procedures. 270 Finally, we outline the stable and time-reversible XL-BOMD algorithm, which ensures energy conservation and dissipates numerical error while simultaneously avoiding the overhead of an iterative SCF optimization required in regular BOMD. 279 6.2.1 Semi-Empirical Quantum Chemistry Methods All semi-empirical methods implemented in PYSEQM use atomic orbitals (AOs, 𝜙 𝜇 (𝒓 ) ) as basis functions, with most elemental basis sets only containing valence shell AOs. Slater-type AOs are used by all of these methods and are detailed here: 𝜙 𝑛𝑙𝑚 =𝑅 𝑛𝑙 (𝑟 )𝑦 𝑙𝑚 (𝘃 ,𝜙 ) (6.1) 𝑅 𝑛𝑙 (𝑟 )=(2𝘁 𝑛𝑙 ) 𝑛 +1/2 [(2𝑛 )!] −1/2 𝑟 𝑛 −1 𝑒 −𝜁 𝑛𝑙 𝑟 (6.2) 83 where n, l and m are quantum numbers for 𝜙 𝜇 (𝒓 ) , 𝑅 𝑛𝑙 (𝑟 ) is the radial function with orbital exponent 𝘁 𝑛𝑙 and 𝑦 𝑙𝑚 (𝘃 ,𝜙 ) is the real spherical harmonic function. The molecular orbitals (MOs, 𝜓 𝑖 (𝒓 ) ) are expressed as a linear combination of AOs: 𝜓 𝑖 (𝑟 )=∑𝐶 𝑖 𝜇 𝜙 𝜇 (𝑟 ) 𝜇 (6.3) The expansion coefficients 𝐶 𝑖𝜇 are obtained by solving the Roothaan-Hall equation: 288 𝑭𝑪 =𝑺𝑪𝑬 (6.4) where 𝑭 is the Fock Matrix, 𝑪 is the matrix of 𝐶 𝑖𝜇 , 𝑬 is the diagonal matrix with MO energies, 𝑺 is the overlap matrix. These semi-empirical methods make the Zero-Differential Overlap (ZDO) approximation, 252 setting 𝑺 =𝑰 reducing the equation to: 𝑭𝑪 =𝑪𝑬 (6.5) 𝑭 is composed with one-electron (𝒉 ) and two-electron (𝑮 ) parts, which for closed shell systems are expressed as: 𝑭 (𝑫 )=𝒉 +𝑮 (𝑫 ) (6.6) ℎ 𝜇𝜈 =⟨𝜇 |− 1 2 𝛻 2 |𝜈 ⟩−∑𝑍 𝐴 ⟨𝜇 | 1 𝑅 𝐴 |𝜈 ⟩ 𝐴 (6.7) 𝐺 𝜇𝜈 =∑𝐷 𝜆𝜎 [(𝜇𝜈 ,𝜆𝜎 )− 1 2 (𝜇𝜆 ,𝜈𝜎 )] 𝜆𝜎 (𝜇𝜈 ,𝜆𝜎 )=𝑒 2 ∬𝜙 𝜇 (𝒓 𝟏 )𝜙 𝜈 (𝒓 𝟏 ) 1 𝑟 12 𝜙 𝜆 (𝒓 𝟐 )𝜙 𝜎 (𝒓 𝟐 )d𝒓 𝟏 d𝒓 𝟐 (6.8) where 𝑍 𝐴 is the effective charge for nuclei A with core shell electrons, 1 𝑅 𝐴 is the one electron potential energy operator, and − 1 2 𝛻 2 is the kinetic energy operator. The density matrix 𝐷 𝜆𝜎 can be computed from the molecular orbital coefficients of the occupied states in the following way: 84 𝐷 𝜆𝜎 =2∑𝐶 𝑖𝜆 𝐶 𝑖𝜎 𝑖 (6.9) The MNDO, AM1 and PM3 methods utilize the Neglect of Differential-Diatomic Overlap (NDDO) approximation, where integrals 𝜙 𝜇 (𝒓 )𝜙 𝜈 (𝒓 )d𝒓 are not vanishing only if 𝜙 𝜇 (𝒓 ) and 𝜙 𝜈 (𝒓 ) are centered on same atom. With NDDO, all three- and four-center integrals in G are ignored, reducing the total number of integrals from O(N 4 ) in G to O(N 2 ), where N is the total number of atoms. With the NDDO approximation in the MNDO-based model, the one-electron matrix 𝒉 is further approximated as: ℎ 𝜇𝜈 = { 𝑈 𝜇𝜇 −∑𝑍 𝐵 (𝜇𝜇 ,𝑠 𝐵 𝑠 𝐵 ) 𝜇 =𝜈 𝐵 ≠𝐴 −∑𝑍 𝐵 (𝜇𝜈 ,𝑠 𝐵 𝑠 𝐵 ) 𝐵 ≠𝐴 𝜇 ,𝜈 centered on atom 𝐴 𝛽 𝜇𝜐 = (𝛽 𝜇 +𝛽 𝜈 ) 2 𝑆 𝜇𝜈 otherwise (6.10) where |𝑠 𝐵 ⟩ is the s-orbital of valence shell for atom B and the two-center one-electron integral from the interaction with core shell of atom B is evaluated as ⟨𝜇 |1 𝑅 𝐵 ⁄ |𝜈 ⟩=(𝜇𝜈 ,𝑠 𝐵 𝑠 𝐵 ) if 𝜇 ,𝜈 are centered on same atom. 𝑈 𝜇𝜇 =⟨𝜇 |−𝛻 2 2 ⁄ |𝜇 ⟩−𝑍 𝐴 ⟨𝜇 |1 𝑅 𝐴 ⁄ |𝜇 ⟩ is the onsite energy for orbital 𝜙 𝜇 on atom A, and is parameterized for each element. The resonance integrals 𝛽 𝜇𝜐 are approximated as the average of each orbital’s β parameter multiplied by the corresponding term in the overlap matrix. With NNDO, the two-electron integral (𝜇𝜈 ,𝜆𝜎 ) in G (Eq. (6.8)) are not be ignored only if both pairs of orbitals (𝜙 𝜇 ,𝜙 𝜈 ) and (𝜙 𝜆 ,𝜙 𝜎 ) are located on the same atoms, respectively. For evaluating (𝜇𝜈 ,𝜆𝜎 ) , a classical approximation is used. (𝜇𝜈 ,𝜆𝜎 ) describes the electron interaction energy between charge distribution 𝑒𝜙 𝜇 𝜙 𝜈 at atom A and 𝑒𝜙 𝜆 𝜙 𝜎 on atom B. e𝜙 𝜇 𝜙 𝜈 is represented 85 by a set of multipoles {𝑀 𝑙𝑚 𝐴 }, and each multipole is represented by a configuration of set of point charges with net charge e and multipole 𝑀 𝑙𝑚 𝐴 . 289 The integrals are evaluated as (𝜇𝜈 ,𝜆𝜎 )=∑∑∑[𝑀 𝑙 1 𝑚 𝐴 ,𝑀 𝑙 2 𝑚 B ] 𝑚 𝑙 2 𝑙 1 = 𝑒 2 2 𝑙 1 +𝑙 2 ∑∑𝑓 (𝑅 𝑖𝑗 ) 2 𝑙 2 𝑗 =1 2 𝑙 1 𝑖 =1 (6.11) where l is the multipole order and m is the multipole tensor index, and 𝑓 (𝑅 ) is Coulomb interaction between each unit point charge. To satisfy the asymptotic behavior with 𝑅 →0 + ,+∞, Dewar- Sabelli-Klopman (DSK) approximation 254, 289 is used: 𝑓 (𝑅 𝑖𝑗 )=[𝑅 𝑖𝑗 2 +(ρ 𝑙 1 𝐴 +ρ 𝑙 2 𝐵 ) 2 ] −1/2 (6.12) where 𝜌 𝑙 𝐴 is the additive term for atom A with multipole order l. The value of 𝜌 𝑙 𝐴 is derived for each atomic species in the limit 𝑅 →0 + for same type of atom A and B, to yield the correct direct and exchange integrals 𝑔 𝜇𝜈 , ℎ 𝜇𝜈 : lim 𝑅 𝐴𝐵 →0 + (𝜇 𝐴 𝜇 𝐴 ,𝜈 𝐵 𝜈 𝐵 )=(𝜇 𝐴 𝜇 𝐴 ,𝜈 𝐴 𝜈 𝐴 )=𝑔 𝜇𝜈 (6.13) lim 𝑅 𝐴𝐵 →0 + (𝜇 𝐴 𝜈 𝐴 ,𝜇 𝐵 𝜈 𝐵 )=(𝜇 𝐴 𝜈 𝐴 ,𝜇 𝐴 𝜈 𝐴 )=ℎ 𝜇𝜈 (6.14) here 𝑔 𝜇𝜈 , ℎ 𝜇𝜈 are parameterized from experimental or ab initio calculation data, and they are smaller than the analytic values as to compensate for the neglect of there- and four-center coulomb integrals and the exclusion of electron correlation in HF. 254-255 Apart from the parameterization in the Hamiltonian, the only difference between MNDO, AM1, and PM3 is the nuclear-nuclear interaction term. The nuclear energy between atom A and B in MNDO is 𝐸 𝑛𝑢𝑐 𝐴𝐵 =𝑍 A 𝑍 B (𝑠 A s A ,s B s B )[1+𝑒 −𝛼 A 𝑅 AB +𝑒 −𝛼 B 𝑅 AB ] (6.15) 86 where 𝛼 A ,𝛼 B are atomic parameters for atom A and B respectively. 𝑅 AB is simply the distance between the two atoms. In contrast, AM1 and PM3 have the following nuclear-nuclear interaction term: 𝐸 𝑛𝑢𝑐 𝐴𝐵 =𝑍 𝐴 𝑍 𝐵 (𝑠 𝐴 𝑠 𝐴 ,𝑠 𝐵 𝑠 𝐵 )[𝑒 −α 𝐴 𝑅 𝐴𝐵 +𝑒 −α 𝐵 𝑅 𝐴𝐵 +𝐹 𝐴 (𝑅 𝐴 𝐵 )+𝐹 𝐵 (𝑅 𝐴𝐵 )] (6.16) 𝐹 𝐴 (𝑅 𝐴𝐵 )=∑𝐾 𝐴 𝑖 exp[𝐿 𝐴 𝑖 (𝑅 𝐴𝐵 −𝑀 𝐴 𝑖 ) 2 ] 𝑚 𝐴 𝑖 (6.17) where 𝐾 𝐴 𝑖 , 𝐿 𝐴 𝑖 and 𝑀 𝐴 𝑖 are sequences of atom-type dependent parameters for the Gaussian expansion of 𝐹 𝐴 (𝑅 𝐴𝐵 ) . In AM1 there are 2 to 4 terms in this sum, depending on the element, while for PM3 there are only 2. With this, the total energy of the system is then 𝐸 𝑡𝑜𝑡 =𝐸 elec +∑𝐸 𝑛𝑢𝑐 𝐴𝐵 𝐴 <𝐵 (6.18) where the electronic energy is 𝐸 𝑒𝑙 ec (𝑹 ,𝑫 )= 1 2 Tr[𝑫 ∗(𝒉 +𝑭 (𝑫 ))] (6.19) Eq. (6.5) is solved iteratively using an SCF loop. A trial density matrix 𝑫 0 is constructed with the valence charge of each atom evenly distributed on its valence shell orbitals. Other strategies of initialization 𝑫 0 adopted in other packages includes using a faster QM method first, like the Hückel model, or evenly distribution of the electron denisty with a small noise. During each step k, 𝑫 𝑘 −1 is used to form the Fock matrix 𝑭 (𝑫 𝑘 −1 ) based on Eq. (6.6), and the density 𝑫 𝑘 ̃ is constructed either by diagonalizing F(𝑫 𝑘 −1 ) or by some expansion scheme as discussed below. The new density for the next step 𝑫 𝑘 is formed by linear mixing of 𝑫 𝑘 ̃ ,𝑫 𝑘 −1 ,𝑫 𝑘 −2 ,...𝑫 𝑘 −𝑚 , with the mixing coefficients being evaluated by linear interpolation algorithms like adaptive mixing 290 and Pulay 291 . Alternatively, it can be constructed with algorithms such as Camp-King 292 , or SOSCF. 293 The loop can be stopped with a variety of different convergence criterion, such as a 87 small ∆𝐸 𝑡𝑜𝑡 between subsequent iterations, a small ∆𝑫 between subsequent iterations, or when the density matrix and the Fock matrix commute to within some threshold 𝜺 ([𝑫 𝑘 ,𝑭 (𝑫 𝑘 )]= 𝑫 𝑘 𝑭 (𝑫 𝑘 )−𝑭 (𝑫 𝑘 )𝑫 𝑘 <𝜺 ). The forces acting on the atoms are further computed as the negative gradient of total energy with respect to the nuclear coordinates. For BOMD, the electrons are assumed to stay in the ground electronic state, and due to Hellmann-Feynman theorem the gradient on the density matrix is assumed to be zero. Using PyTorch’s Automatic Differentiation features, 271 the analytic gradient is automatically obtained as the total energy is computed. Automatic Differentiation does not rely on finite differences to compute gradients, instead constructing a computational graph and saving the required intermediate variables for the backward gradient calculations, using the chain rule to calculate gradients in an efficient manner. 6.2.2 SP2 algorithm For each SCF loop, a new density matrix 𝑫 𝑘 ̃ must be constructed from the Fock Matrix 𝑭 𝑘 . Traditionally, this is done by diagonalizing 𝑭 𝑘 and forming 𝑫 𝑘 ̃ using the eigenvectors based on Eq. (6.9). As an alternative approach, we can directly form 𝑫 𝑘 ̃ from 𝑭 𝑘 without matrix diagonalization using the property that 𝑫 𝑘 ̃ is given by a step function of 𝑭 𝑘 with the step formed at the chemical potential, as is described in Eq. (6.21) below. A step function expansion can then be performed using various recursive Fermi-operator schemes that avoid any explicit diagonalization. The simplest and possibly most efficient scheme is based on second-order spectral projection polynomials, known as the SP2 algorithm. The computational kernel of the SP2 scheme that dominates the cost is a generalized matrix-matrix multiplication, which can make efficient use of hardware acceleration, such as GPUs. 270 To begin, 𝑫 𝑘 ̃ can be rewritten as: 𝑫 𝑘 ̃ =𝘃 (𝜇 𝑰 −𝑭 𝑘 ),Tr[𝑫 𝑘 ̃ ]=𝑁 𝑜𝑐𝑐 (6.20) 88 where 𝘃 (𝑥 ) is the Heaviside step function, 𝜇 and Nocc are the chemical potential and the occupation number of the system. Generally, 𝘃 (𝑥 −𝜇 ) can be expanded and approximated by a serial polynomial expansion, e.g., with Chebyshev polynomials, with the step centered on a predefined 𝜇 . 294 The SP2 algorithm approximates the step function with recursive expansion based on second- order spectral projection polynomials given in Eq. (6.24), which have the stationary end points at 0 and 1, without requiring prior knowledge of 𝜇 . In the SP2 scheme the Fock matrix 𝑭 𝑘 is first scaled so the eigen-spectrum lies on the interval [0,1]: 𝑿 (𝟎 ) = (𝘀 𝑁 𝑰 −𝑭 𝑘 ) (𝘀 𝑁 −𝘀 1 ) (6.21) where 𝘀 1 and 𝘀 𝑁 are estimates of the smallest and largest eigenvalues of F, respectively. 𝘀 1 and 𝘀 𝑁 can be estimated, for example, with Gershgorin Circle Theorem 294 , which dictates that all the eigenvalues of the Hermitian matrix F lie in the interval on real axis, [𝐦𝐢𝐧 𝑖 (𝑅 𝑖 −𝑭 𝑖𝑖 ),𝐦𝐚𝐱 𝑖 (𝑅 𝑖 +𝑭 𝑖𝑖 )],𝑅 𝑖 =∑|𝑭 𝑖𝑗 | 𝑗 ≠𝑖 (6.22) where 𝑅 𝑖 and 𝑭 𝑖𝑖 are the radius and center of Gershgorin Disc i. The estimated bounds are used to replace 𝘀 1 and 𝘀 𝑁 in the scaling of F as 𝘀 1/𝑁 are not known a priori. A quadratic form of the spectral projection (or purification) polynomials P (a) and P (b) are defined as: 𝑃 (𝑎 ) (𝑋 )=𝑋 2 ,𝑃 (𝑏 ) (𝑋 )=2𝑋 −𝑋 2 (6.23) These projection polynomials are applied iteratively to 𝑿 (𝟎 ) to generate a converging sequence of 𝑿 (𝒏 ) under the following rules: 𝑿 (𝒏 +𝟏 ) ={ 𝑃 (𝑎 ) (𝑿 (𝒏 ) ),|Tr[𝑃 (𝑎 ) (𝑿 (𝒏 ) )]−𝑁 𝑜𝑐𝑐 |< |Tr[𝑃 (𝑏 ) (𝑿 (𝒏 ) )]−𝑁 𝑜𝑐𝑐 | 𝑃 (𝑏 ) (𝑿 (𝒏 ) ),otherwise (6.24) 89 In this way, the eigenvalues of 𝑿 (𝒏 ) converge to the step function of H as 𝑛 →∞, while Tr(𝑿 (𝒏 ) ) converges to the occupied number of orbitals 𝑁 𝑜𝑐𝑐 . The iterative procedure is stopped when the error in total electron count at step n and n-1 is smaller than a 10 -4 . The error in electron number at step n is defined as: err(𝑛 )=|Tr(𝑿 (𝒏 ) )−𝑁 𝑜𝑐𝑐 | (6.25) which is used as a convergence test together with the idempotency error that can be estimated by Tr[P(I-P)]. Typically, 20 to 30 iterations are required before a sufficient convergence is achieved. Finally, the density matrix can be formed from the converged 𝑿 (𝒏 ) 𝑫 𝑘 ̃ =2𝐥𝐢𝐦 𝑛 →∞ 𝑿 (𝒏 ) (6.26) where the factor of 2 accounts for the spin degeneracy in a closed shell system. 6.2.3 Extended Lagrangian BOMD In conventional Born-Oppenheimer molecular dynamics, the motion of the nuclei is described by the Lagrangian: ℒ 𝐵𝑂 (𝑹 ,𝑹 ̇)= 1 2 ∑𝑚 𝑘 𝑅 ̇ 𝑘 2 −𝑈 (𝑹 ,𝑫 ) 𝑘 (6.27) Here 𝑈 (𝑹 ,𝑫 ) is the potential energy surface with ground state density matrix 𝑫 . For a given 𝑹 , U is the total energy defined in Eq. (6.19). In regular BOMD, an extrapolation of ground state density matrices from previous time steps are used as initial guess to the SCF optimization procedure. This extrapolation, followed by an SCF optimization, which in practice is never complete, breaks the time-reversibility of the underlying propagation of the electronic degrees of freedom. This common form of direct BOMD simulation typically leads to an unphysical systematic drift in the total energy and phase space. 295 The problem can be reduced or removed by tightening the SCF convergence or by ignoring the extrapolations and instead starting each new SCF optimization 90 from scratch using, for example, overlapping atomic densities in each new time step. This may substantially increase the computational cost. To avoid these shortcomings of regular direct BOMD, a time reversible formulation can be constructed by introducing an auxiliary electronic density matrix 𝑷 as a dynamical tensor variable that evolves in addition to the nuclear coordinates and their velocities. This can be achieved using XL-BOMD 278-282 , which is defined through an extended Lagrangian, ℒ 𝑋𝐿𝐵𝑂 (𝑹 ,𝑹 ̇,𝑷 ,𝑷 ̇)= 1 2 ∑𝑚 𝑘 𝑅 ̇ 𝑘 2 −𝑈 (𝑹 ,𝑷 ) 𝑘 + 𝜇 2 Tr[𝑷 𝟐̇ ]− 𝜇 𝜔 2 2 Tr[(𝑫 −𝑷 ) 𝑇 𝑲 𝑇 𝑲 (𝑫 −𝑷 )](6.28) Here 𝜇 is a fictitious electron mass parameter, 𝜔 is the frequency parameters for the electronic degree of freedom and 𝑲 𝑇 𝑲 is a metric tensor of the harmonic well. These variables dictate how the auxiliary electronic density matrix 𝑷 , fluctuates through a harmonic oscillator located around an optimized electronic density matrix 𝑫 . The modified “shadow” potential energy surface 𝑈 (𝑹 ,𝑷 ) is defined through a constrained density matrix minimization as: 𝑈 (𝑹 ,𝑷 )=𝐸 (1) 𝑒𝑙𝑒𝑐 (𝑹 ,𝑷 )+𝐸 𝑛𝑢𝑐 (𝑹 ) (6.29) 𝐸 (1) 𝑒𝑙𝑒𝑐 (𝑹 ,𝑷 )= 1 2 (2Tr[𝒉𝑫 (𝑷 )]+𝑇𝑟 [(𝟐𝑫 (𝑷 )− 𝑷 )𝑮 (𝑷 )]) (6.30) The P dependency of D(P) is dropped in most part of the text for simplicity. It is shown by Niklasson and Cawkwell, 281 the difference between 𝑫 and 𝑷 scales as (𝑫 −𝑷 )= 𝓞 (Ω 2 /𝜔 2 ) , where Ω is the highest characteristic frequency of nuclear motion. In an adiabatic limit where 𝜇 → 0, with 𝜇𝜔 being constant, assuming 𝜔 ≫Ω, and with K chosen as a super operator given by the inverse Jacobian 281 of the residual D-P, then the Euler-Lagrange equations of motion for ℒ 𝑋𝐿𝐵𝑂 that governs the time evolution of 𝑷 and 𝑹 , is given by 𝑚 𝑘 𝑅 𝑘 ̈ =− ∂𝑈 (𝑹 ,𝑷 ) ∂𝑅 𝑘 = 1 2 Tr[2 𝜕 𝒉 ∂𝑅 𝑘 𝑫 +(𝟐𝑫 −𝑷 ) 𝜕 𝑮 (𝑷 ) ∂𝑅 𝑘 ] (6.31) 91 𝑷 ̈ =−ω 2 𝑲 (𝑫 −𝑷 ) (6.32) where 𝑫 (𝑡 ) is obtained at the constrained minimization, as the lowest stationary solution of the linearized energy expression in Eq. (6.31). Thanks to the linearization, the minimization is achieved in a single step and no iterative SCF optimization is involved. In all of our calculations we will use a scaled delta function approximation of the kernel K, where is replaced by -c*I, with c in the interval [0,1]. The initial value of 𝑷 (𝑡 ) is given by a full regular SCF optimization of 𝑫 (0) , i.e. where 𝑷 (0) = 𝑫 (0) with 𝑷 ̇(0) =𝟎 . The nuclear equations of motions are integrated with a symplectic velocity Verlet scheme. 296 Eq.(6.33) shows that 𝑷 oscillates on a harmonic surface centered around the optimized density matrix 𝑫 . In the adiabatic limit(𝑫 −𝑷 )= 𝓞 (Ω 2 /𝜔 2 ). The equations of motion for the electronic degrees of freedom, 𝑷 , can be integrated using the standard Verlet scheme, 296 𝑷 (𝑡 +∆𝑡 )=2𝑷 (𝑡 )−𝑷 (𝑡 −∆𝑡 )+∆𝑡 2 𝑷 ̈. (6.33) As the integration is time reversible, the SCF optimization in Eq. (6.35) keeps the time-reversal symmetry in the evolution of 𝑫 (𝒕 ) . The systematic drift of total energy and in phase space due to breaking time-reversal symmetry in BOMD can thus be eliminated with the XL-BOMD framework. However, the time-reversible integration of XL-BOMD introduces a problem due to internal numerical noise from an approximate matrix algebra based on finite floating-point operations. 279 A time-reversible integration scheme does not dissipate any numerical noise. Numerical error therefore accumulates and may eventually lead to a significant deviation of P from D. This subsequently results in a drift of the potential energy surface 𝑈 (𝑹 ,𝑷 ) when performing MD simulations for long periods of time. To avoid this noise accumulation, a dampening force term is added to the Verlet integration of 𝑷 (𝑡 ) : 92 𝑷 (𝑡 +∆𝑡 )=2𝑷 (𝑡 )−𝑷 (𝑡 −∆𝑡 )+∆𝑡 2 𝑷 ̈ +α∑𝑐 𝑘 𝑷 (𝑡 −𝑘 ∆𝑡 ) 𝐾 𝑘 =0 . (6.34) The dissipation is introduced by the last summation term that dampens out noise accumulation. Here 𝛼 , and 𝑐 𝑘 are parameters with their values chosen such that the error is reduced while keeping the time-reversibility to a higher odd order, O(∆𝑡 (2𝐾 −3) ), of the integration time step. Their optimized values can be found in Ref [ 279 ] together with the value of 𝜅 =∆𝑡 2 𝜔 2 . The integration scheme dissipates numerical errors and generates stable MD trajectories with a well-controlled long-term energy conservation. 6.3 Results and Discussion To test the models and algorithms implemented in the PYSEQM package and benchmark performance, we run conventional BOMD and XL-BOMD on two model systems. The first system, a “nanostar” phenylene-ethynylene dendrimer with 884 atoms, 287 is used as a showcase molecule. The second system, 8 polyethylene molecules with sizes ranging from H-(C2H4)1-H to H- (C2H4)128-H, is used to illustrate the scaling of the code with a system size. Benchmarks are performed on both CPU and GPU architectures with double precision. To put these benchmarks in perspective of conventional implementation of SEQM, the results are compared with a standard quantum chemistry package (ORCA). 297 As PyTorch relies on an asynchronous execution when utilizing GPU, synchronizing calls are used to enforce accurate GPU timings, and based on the performance, adding synchronizing slows down the performance by around 2%. 6.3.1 Nanostar Dendrimer The nanostar tested here is a dendrimer with two unique chemical environments, a core ethynylperylene chromophore and phenylene-ethynylene terminals. 287 Because of potential wide applications and fruitful experimental and theoretical studies on this type of dendrimer, 287, 298-299 93 the nanostar, containing 460 C and 424 H atoms in total, makes a good benchmark system due to its size. The structure was initially equilibrated at 300K for 50 ps using the AIREBO force field 300 in the LAMMPS molecular dynamics (MD) package. 48 Initializing from the equilibrated structure, and randomly assigned velocities, we further run MD for 200 steps with a 0.2 fs timestep. Simulations for all combinations of BOMD/XL-BOMD, SP2/conventional diagonalization, and CPU/GPU are performed to get a comprehensive measure of code performance under a variety of different simulation conditions. The simulations are tested on one Intel Xeon E5-2660_v3 CPU core and one Nvidia Tesla V100-SXM2 GPU. The convergence criteria for SCF procedure in BOMD is set to be 10 -6 Eh (Hartree Energy). Finally, a BOMD simulation performed with ORCA is run for a timing comparison. Figure 6.1 shows the average time cost per MD step in both BOMD and XL-BOMD. Additionally, these timings are reported separately for CPU and GPU execution. Modern GPUs are capable of operating in “single instruction multiple data” (SIMD) mode across thousands of cores on single GPU, resulting in a factor of 10 to 20 speed up when operating on a GPU vs. a CPU, as shown in Fig. 6.1(a) and (b). In Fig. 6.1(a), we report the execution time of a BOMD simulation performed on a CPU in both PYSEQM and ORCA software. The PYSEQM calculation is approximately 1.5 times faster than ORCA even though ORCA utilizes a more sophisticated SCF algorithm that generally requires fewer iterations to achieve convergence than the SCF algorithm used in PYSEQM. The improvement is mainly coming from the SCF procedure, as computing force, Hcore and Coulombic Integrals takes approximately same amount of time between the two codes. This is due to implementation details: ORCA utilizes the Second Order SCF orbital optimization (SOSCF) 301 algorithm, which is more numerically stable (see our discussion in the next section) and costly. The time spent in the SCF cycle is greatly reduced by 94 the XL-BOMD algorithm because it only requires a single matrix diagonalization to obtain the correct density. This leads to approximately a 95% reduction in compute time for CPU calculations and an 85% reduction for GPU computations for the SCF loop. This translates to an overall factor of 2 speedup for the full MD propagation. The SP2 algorithm performs differently on CPU and GPU: it is slightly slower on the CPU, but 3 times faster on the GPU. The SP2 algorithm performs much better on GPU architectures because the dominant operations in this algorithm are matrix- matrix multiplications, which benefit greatly from GPU acceleration. Overall, PYSEQM utilizing the XL-BOMD and SP2 algorithms running on a GPU is more than 50 times faster than ORCA running on a CPU for the nanostar system. Figure 6.1: Stacked Histogram of the average time spent computing forces (red), Hcore and integrals (green), the SCF loop (blue), and other code (white) per MD step with the nanostar dendrimer. This is both with and without SP2 and the XL method on an Intel Xeon E5 CPU (a) and a Nvidia Tesla V100 (b). The left most column is a BOMD computed with ORCA on the CPU for comparison. It is critical that quantum-based MD simulations conserve total energy when simulating an NVE ensemble. In practice, most of such simulations experience some drift in the total energy over time, due in part to the finite convergence of the SCF procedure. To verify that PYSEQM is capable of conserving energy to within normal tolerances, we performed a 0.6 ps NVE MD 95 simulation on the nanostar system with a variety of algorithms and both a 0.1 and 0.2 fs time step. The total energy of these simulations is shown in Fig. 6.2. Table 6.1 shows the drifting coefficient, which is an average per atom unit of energy change per unit time, and energy fluctuation, which is the standard deviation of total energy after removing the drift. When restarting the SCF procedure with the density matrix obtained from the previous MD step, a large drift in total energy is observed as seen in the yellow trace of Fig. 6.2. This can be improved in BOMD simulations with sophisticated SCF algorithms like SOSCF, 301 Direct Inversion in the Iterative Subspace (DIIS), 302 which lead to much better energy conservation, for example, in ORCA. In PYSEQM, the drift can also be fixed through the use of the XL-BOMD algorithm, which for this particular system, has a drift of 2.94  10 -3 meV/ps/atom with dt = 0.1 fs. This amount of drift is comparable to a classical MD simulation performed with LAMMPS using the AIREBO (not shown), 48, 300 which has a drift around 1.0  10 -3 meV/ps/atom with dt = 0.1 fs. The energy drift in PYSEQM BOMD can also be reduced by restarting the SCF procedure from a new guess at each time step (Fig. 6.2, red right triangle trace) rather than re-using the previous time step’s converged density matrix, though this is not used in practice for two reasons. First, starting from a guess density at each time step greatly increases the number of SCF iterations required at each point in the MD trajectory. Second, the density matrix may become stuck in local minima during the SCF procedure, resulting in a small energy jump as shown in Fig. 6.2 (red right triangle trace). The flat plateau of the curves in Fig. 6.2 (red right triangle trace) shows that energy conservation in BOMD simulations can be improved by not reusing the density matrix when compared to the results in Fig. 6.2 (yellow star trace). In practice, a general acceptable drift is around 1 meV/ps/atom over long time scales, and the discussed strategy is not adapted for efficiency. 96 Another issue with respect to the energy conservation is the fluctuation of energy. This is coming from the discretization of propagation. Under the Velocity Verlet scheme, the fluctuation is in order of O(∆𝑡 2 ). 303-304 Because of this, the fluctuation can be suppressed by reducing the timestep used for propagation. However, in practice for efficiency, the timestep is frequently chosen to be 0.05 to 0.01 times the timescale for the highest characteristic vibrational frequency of the system. Generally, this guideline results in a timestep of 0.1 fs to 1 fs for organic systems. Here with reducing timestep from 0.2 fs to 0.1 fs, the energy fluctuation is reduced by around 50% with BOMD and more than 75% with XL-BOMD. Figure 6.2: Total energy drift of the nanostar system during 0.6 ps simulations under various simulation conditions, 0.2 ps is shown here. The time steps (a) dt = 0.2 fs and (b) 0.1 fs are used with BOMD (yellow star), BOMD without reusing density matrix in consecutive MD steps (read right triangle), XL-BOMD (blue left triangle), BOMD on ORCA (black dot). 97 BOMD BOMD XL-BOMD XL-BOMD BOMD(ORCA) Timestep (fs) 0.2 0.1 0.2 0.1 0.2 Drift (meV/ps/atom) -2.13 -2.03 -1.21e-2 -2.94e-3 -1.95e-1 Fluctuation (eV) 1.63e-2 1.17e-2 1.77e-2 4.42e-3 1.49e-2 Table 6.1: Energy drift and fluctuation from Fig. 6.2 for BOMD and XL-BOMD done with package PYSEQM and ORCA. BOMD results here are for the cases with re-using density matrix in consecutive MD steps. 6.3.2 Polyethylene As the simplest polymer, polyethylene is a good system with which to check the scaling of PYSEQM with molecular size. We create 8 polyethylene chains H-(C2H4)n-H with n = 1,2, 4,…, 128, with the largest having a total 770 atoms (256 C and 512 H). Similar to the nanostar, these molecules are first relaxed and equilibrated at 300K with LAMMPS using the AIREBO force field for 50 ps. Then we perform simulations for 20 steps with a timestep of 0.2 fs, with all combinations of BOMD/XL-BOMD, SP2/conventional diagonalization, and CPU (Intel Xeon E5)/GPU (Nvidia TITAN_V) for these 8 molecules. This gives a good benchmark to quantify the scaling of PYSEQM with system size. Figure 6.3 shows the scaling of the Hcore and Coulombic integrals, Force, SCF and MD computations. When running on a CPU, ORCA is faster for small molecules, but it scales worse than PYSEQM as shown in Fig. 6.3(a)-(d), with ORCA’s performance becoming equal to that of PYSEQM for molecules n>=32 as shown in Fig. 6.3(d). It is reasonable that PYSEQM is slower for smaller systems because the overhead for function calls dominate the time cost in PYSEQM for small systems. This overhead in python is the main reason why it is slower than other low- level programming language such as C/C++ used by ORCA. However, for larger systems, where 98 overhead takes a smaller fraction of the total time of simulation, PYSEQM outperforms ORCA because it uses a vectorization strategy for computing Hcore, Integrals, and Force as shown in Fig. 6.3(a) and (b). Figures 6.3(a) and (b) show results from using PYSEQM on CPU and GPU and ORCA on CPU. For small molecules, ORCA performs worse than PYSEQM in the SCF procedure, which is because ORCA uses a systematic combination of SCF algorithms: in general, it is more stable and robust but slower. The difference with using SP2 or conventional diagonalization on CPU is negligible and the results are not shown in Fig. 6.3. Combined with the results on the nanostar molecule, PYSEQM is faster than ORCA for larger molecules when using a CPU. However, the true strength of PYSEQM is its ability to run on a GPU. This results in a significant gain in performance when running on large systems. For small systems, the GPU mode is less efficient. Modern GPUs have thousands of slow cores, which only perform well with large amounts of data being processed by the same operations. As a result, for small systems the GPU hardware is not fully utilized. In turn, the GPU mode of PYSEQM shows poor scaling for small systems, with the total time required to perform a calculation nearly independent of systems size for small molecules, as shown in Fig. 6.3. However, the ability of PYSEQM to run multiple simulations in a batch mode (discussed below) compensates for this drawback. Once the system size is large enough to use a significant fraction of the GPU, PYSEQM on a GPU can easily outperform PYSEQM and ORCA on CPU. There is similar trend for using the SP2 algorithm on a GPU. For smaller systems, SP2 is slower than conventional diagonalization, but for large systems, the matrix multiplications required by SP2 can fully utilize GPU cores and outperform conventional diagonalization algorithms. 99 Figure 6.3: Average time spending per MD step on computing (a) Hcore and Coulombic Integrals, (b) Force, (c) SCF procedure and (d) MD for eight molecules H-(C2H4)n-H with size n = 1,2, …, 128. The timing is done on GPU (blue), CPU (red) and CPU with ORCA (black dot) with BOMD (no marker) and XL-BOMD (vertical marker |) with SP2 (solid line) or conventional diagonalization (dash line). Only the results with SP2 on CPU is shown here for clearness of figures as the difference with or without using SP2 is small. Due to our desire to interface PYSEQM with machine learning algorithms, it has the ability to run in “batched mode,” where simulations are performed simultaneously on an ensemble of molecules to fully take advantage of the parallel GPU architecture. In addition to being useful for ML applications, many simulations in physics and chemistry require ensemble averaging such as computing the decay time of photoexcited electrons, 305 or estimating quantum yields, 306 making this “batched mode” directly relevant for applications as well. Thus in Fig. 6.4 we benchmark PYSEQM in the batched mode. We run the MD simulations with 32 molecules where trajectories 100 were initiated with different initial velocities and starting from the same nuclei coordinates for each of 8 polyethylene systems. The simulations here are done with the same setting as for the single molecule mode. Figure 6.4 shows the time cost per MD step per molecule when running with batches while Fig. 6.5 shows the relative performance boost of batch mode versus running the simulations in serial. Even running MD with 32 molecules at the same time, the occupancy on GPU is still relatively low for small systems and the scaling is poor as shown in Fig. 6.4. However, the time cost per molecule is greatly reduced not only with GPU, but also when running on CPU. For small systems, the performance of PYSEQM in batched mode on CPU is 10~20 times better than that in the serial mode. As system size grows, the performance boost is quickly degraded on the CPU, because PyTorch switches to perform operations in series for large operations such as matrix diagonalizations or multiplications. This can be seen in the scaling of computing forces and the SCF procedure in Fig. 6.4 (b) and (c). The differences with SP2 and conventional diagonalization are relatively small and thus are not shown in Figs 6.4 and 6.5. On GPU, the performance boost with batch mode is higher and degrades much slower as shown in Fig. 6.5. Additionally, it is around 20~30 times faster than conventional serial simulations on CPU. Additionally, the poor scaling for small systems indicates the further performance gains can be achieved by increasing the batch size to further saturate the GPU. SP2 shows its advantages for n>=8 for our polyethylene systems compared to conventional diagonalization. In our testing with random symmetric matrices, SP2 is faster than direct diagonalization almost for all sizes of matrix. However, as PYSEQM is implemented to be able to deal with ensembles of systems simultaneously, preprocessing is required to use SP2, in which the overhead slows down SP2 for small systems. For larger systems, SP2 is around 5 times faster when compared with using direct 101 diagonalization, as shown in Fig. 6.4 (c), resulting in 2 to 3 times faster in performing MD in Fig. 6.4(d). Overall, from results on Fig. 6.4 (a)-(c), PYSEQM on GPU with batch mode is systematically faster than conventional ORCA calculations by about 20 times. Figure 6.4: Average time needed per MD step per molecule with the batch mode enabled when computing (a) Hcore and Coulombic Integrals, (b) Force, (c) SCF procedure and (d) MD for each of six H-(C2H4)n-H with size n = 1,2, …, 32. Same labeling and coloring are used here as in Fig. 6.3. 102 Figure 6.5: Performance boost when running 32 molecules simultaneously in the batch mode compared with running in serial. Simulations are performed on polyethylene H-(C2H4)n-H, n = 1, 2, …, 32. On CPU (red), the difference is small, and the results with conventional diagonalization are omitted. When computing on GPU (blue), speedups of a factor of 30 are possible for small molecules. 6.4 Conclusions In this work we present the PYSEQM software package, which provides Born Oppenheimer Molecular Dynamics with semi-empirical quantum mechanics methods (MNDO, AM1 and PM3) in the PyTorch framework, enabling GPU acceleration and a native interface with machine learning codes. To ensure rapid ground state SCF computations at each MD steps, PYSEQM implements the SP2 density matrix expansion algorithm. To avoid energy drift present in conventional BOMD based on the quantum methods, PYSEQM further implements XL-BOMD method, which both accelerates the calculations and enforces energy conservation, since an iterative SCF optimization that typically is only approximate is avoided prior to the force evaluations. In particular, this allows achieving the energy conservation 100 to 1000 times better than traditional semi-empirical quantum mechanical molecular dynamics and at a significantly 103 reduction in the computational cost. Finally, to make full use of the highly parallel GPU architecture, PYSEQM supports a “batched mode” where many molecules can be computed simultaneously. This batch feature, where many independent MD trajectories are run simultaneously, is generally absent in many packages, and it is very beneficial for simulations requiring statistical averaging. Overall, the current PYSEQM implementation can treat systems with more than 1000 atoms and with elements from H to Cl, which covers all organic molecules. We further use the nanostar dendrimer and several polyethylene chains to document performance of the PYSEQM code. Various PYSEQM applications to these molecular systems are compared to ORCA computational package featuring conventional implementation to show for the end user PYSEQM’s performance advantages and disadvantages. When running PYSEQM on CPU, its performance is on par with a conventional implementation for smaller systems and is faster for large molecules owing the use of SP2 and XL-BOMD algorithms. Subsequently, here the PYSEQM package provides a valuable alternative to other conventional codes on CPU capable of semiempirical BOMD simulations. In contrast, the main focus of the PYSEQM design is its ability to efficiently utilize GPU. We show that while it is slower to use PYSEQM for small systems on GPU, this can be overcome utilizing the “batch mode”. For large molecular systems (over 500 atoms), PYSEQM becomes very efficient on GPU in serial mode as well. Altogether, we demonstrate a factor of 50 speed up when performing BOMD calculations on GPU when compared to conventional CPU calculations. In summary, in the future we envision multiple applications for and further development of the PYSEQM package. First, this software may be useful for an end user to perform effective BOMD simulations on GPU using tried-and-true old semiempirical models (AM1, MNDO, PM3, and PM6) particularly for large systems. Secondly, our future work will include extensions to other 104 semiempirical methods (overviewed in the Introduction), which will both increase the accuracy as well as facilitate the treatment of more systems. 257 One specific target is extension of the implementation to non-adiabatic excited state molecular dynamics 10, 25, 277 which will enable rapid ensemble propagation using the “batched mode” for calculations such as fewest switches surface hopping 10 (FSSH) or Multiconfigurational Ehrenfest with Ab Initio Multiple Cloning 307 (MCE- AIMC). Finally and most importantly, utilization of the PyTorch framework facilitates the construction of interfaces to machine learning methods, such as HIPNN, 274 ANI, 16, 308 or SchNet 309- 311 . By training such a neural network to generate custom parameter sets for the semi-empirical Hamiltonians, the accuracy of these methods can be greatly increased. This will allow on demand generation of reduced quantum-mechanical models for targeted molecular families able to accurately describe desired properties and dynamics of both ground and excited states. Of particular interest for this method are Hamiltonians utilizing d-oribitals, where we think significant improvements in accuracy can be obtained. 257-258, 263 All these extensions will lead to a wide application of the implementation presented here and be beneficial to researchers working in theoretical computation and machine learning in chemistry, physics and material science. The described PYSEQM package is now released as an open-source software and we invite the community to participate in its further development. 276 105 Chapter 7: Machine Learning with Domain Knowledge in Quantum Chemistry Authors: Guoqing Zhou, Ben Nebgen, Nicholas Lubbers, Sergei Tretiak Contributions: B.N. and S.T. designed the project, G.Z. implemented the model and performed the training and testing and analyzed the results. All contributed to the discussion. 7.1 Introduction Modeling the interacting electrons and nuclei plays a central role in studying chemical and material systems. Traditionally it is done with quantum-mechanical (QM) based level theory like Density Functional Theory (DFT), Coupled Clusters (CC), Configuration Interaction (CI). 30, 167 These methods can provide accurate descriptions of system. However, the poor scaling of their complexity limits the application of these methods in large systems. Apart from that, the material discovery requires to explore the vast chemical compound and configuration space with high efficiency and accuracy. QM-based methods are no longer able to fulfill this gap. The rising of Machine Learning has brought innovation in many areas including material science and chemistry. Two types of strategies have been used to incorporate ML. The first one extract features from the geometries explicitly or implicitly through different structures of neural network, like HIPNN, 274 SchNet, 309 ANI-1. 308 This has been used for predicting various properties like potential energy surface, force, charge, dipoles, etc. And the second one, called Δ-Machine Learning, along with the features derived from geometries, adapts features from the output of low-cost quantum chemistry models, like OrbNet, 312 which use a symmetry adapted atomic orbitals features and lead to a high learning efficiency and extensibility, i.e. applying on different and larger chemical systems. As a comparison, the first strategy requires chunks of data to train, which is usually 106 generated through high-cost DFT methods. And the second one, requiring less data to feed, takes the domain knowledge, and is promising and attracting. However, enormous efforts have been working on the predicting the potential energy surface and force or related properties, 275, 308-309, 312 and lack the flexibility like in the QM-based methods, which provide most desired properties in one calculation. And most critically, these methods show poor extensibility when applying on large systems which are outside the training dataset. As such, we provide a new type of model, and by taking advantage of domain knowledge in quantum chemistry, we build the low-cost Semi-Empirical Quantum Chemistry (SEQM) on top of Hierarchical Interacting Particle Neural Network (HIPNN). 254 The similar strategy has been used with Density Functional Tight Binding (DFTB) and has been proved to be efficient to train with small dataset and increase the accuracy compared to the normal model. 313 Instead of updating charge during the training as in DFTB, we have enabled backpropagation through the Self- Consistent Field (SCF) procedure, and thus make it possible to do the multi-objective training with potential energy surface, force, orbital energies, charge distribution. We train and test models from original PM3 theory, pure HIPNN, and the hybrid model HIPNN+SEQM here on datasets with small organic molecules. And to exam and compare their extensibility, we apply the models on COMP6 dataset, 314 which contains much larger molecules. The results indicate the hybrid model has overall better performance, and at the same time it requires fewer data points for training. We detail the hybrid model structure and the training in section 2. In section 3, we investigate the output of hybrid model, exam and compare these three models’ performance and accuracy. Finally, we summary and conclude in last section and discuss the possible applications. 107 7.2 Model We briefly summarize the structure of HIPNN here as shown in Figure 7.1. A complete description is reported in ref 274 . HIPNN takes molecular configurations as input. Each molecule is represented with a set of atom types and pair distances between atoms, which is invariant under the symmetrical transformations like translation, rotation, reflection, and permutation of atoms. Red blocks in Fig. 7.1 are the on-site layers which apply on the local features for each individual atom. And the green blocks are the interaction layers are applied on the features shared with nearby atoms to transmit information. A sequence of on-site layers and interaction layers are used to process and share features. A fully connected network is then applied on each last on-site layer in the sequences, to yield the deviation from the constant PM3 parameters for the SEQM layers. Figure 7.1: Model structure scheme. Left: HIPNN structure with molecule configuration as input (black and while blocks) and generating features for interacting layers (green blocks) and on-site layers (red blocks) and yield SEQM parameters P with fully connected layers (blue blocks). Middle: SEQM module takes molecule configurations and parameters from HIPNN and generates overlap and Coulomb integrals and performs SCF loop procedure to get various molecular quantities. The SEQM module is detailed in our previous PYSEQM paper. 315 The total energy is expressed as: 108 𝐸 tot =𝐸 elec +∑ 𝐸 nuc 𝐴𝐵 𝐴 <𝐵 (7.1) Where the electronic energy is Eelec = Tr[ D(h + F)]/2, and 𝐸 nuc 𝐴𝐵 is the nuclear interaction between atom A and B. Various empirical parameters are used to construct the h and G and 𝐸 nuc 𝐴𝐵 . In total, for each element in the first three rows of the periodic table, there are 11 parameters for the electronic energy (4 for Hydrogen), and 1 to 13 parameters for the nuclear energy depending on the Semi-Empirical methods (7 for PM3 we used in this article). When interfacing with HIPNN, these semi-empirical parameters are replaced with the dynamically generated ones. We apply softplus (log(1+e x ), here x is the parameter) on certain parameters, if they are required to be positive based on their physics meaning or numerical constraints. The density matrix D is obtained through self-consistent filed (SCF) procedure. SCF requires D used to build the Hamiltonian F should be equal to the one constructed from the occupied eigenvectors of F. This is usually achieved with specific algorithms like adaptive mixing, Pulay, DIIS. 315 Starting with a guess density, D is obtained iteratively. After each step, a new trail D ′ is constructed from D 0 , D 1 , D 2 ,… in previous steps. SCF is stopped until the difference is negatable between two consecutive density matrices. The SCF procedure leads to three problems when interfacing with the neural network part. First, SCF will slow down the training not only on the forward pass, but also on the back- propagating pass. The Hellmann-Feynman theorem allows the gradient computation from total energy(or functions of total energy) to skip the SCF procedure, while this greatly limit the application of SEQM, as many properties are not a function of total energy, like the force which is required for simulating dynamics, or properties like dipoles, which depends on charge distribution. As such, we enable backpropagation through SCF when the target is not limited to total energy related terms. Second, there are many intermediate matrices required to be stored for 109 constructing the computation graph for the backpropagation of SCF. Due to the limitation of device memory, the batch size is an order smaller than the normal neural network model. For these two issues, we make a tradeoff between the learning efficiency and time consumption. A small fraction of molecules is randomly sample from the dataset for training with a small batch size to ensure enough gradient updates. And at the same time, we enforce an upper bound on the number of SCF steps, to reduce the memory usage. The last and most serious problem is that the SCF may fail to converge, or requires too many steps to converge for certain molecules or with certain parameters. The SCF may fail to converge for a large portion of molecules or even all of them in the batch during training. This will make the training unstable and leads to the loss jiggling during the training and make it impossible to update the network parameters. It turns out that the SCF failure happens when the output SEQM parameters from HIPNN may go wildly, far away from the constant PM3 parameter set. For this, we apply a soft constraint on the deviation of SEQM parameters from the HIPNN layers, and at the same time filter out the molecules whose SCF fails with the pre-defined maximal steps. To reduce the jiggling of loss, we scale the loss/gradient by √𝑝 , here p is the fraction of molecules in the batch whose SCF succeed to converge. We scale the loss by √𝑝 assuming the variance for each batch is inversely proportional to number of samples. To accelerate the training, we start by applying a relatively large converging threshold 10 - 5 Eh (Hartree) on the SCF procedure, along the updating of the learning rate, we decay the threshold with factor 0.98 to make it tighter until a minimal value 10 -6 Eh is reached. To constrain the SEQM parameters and minimize the effect from SCF failure, we add the deviation of SEQM parameters to the loss. The total loss is defined as: ℒ=𝑎𝑓 (𝑦̂−𝑦 𝐷𝐹𝑇 )+𝑏 ℒ 𝐿 2 +𝑐𝑔 (∆𝑷 ) (7.2) 110 Here 𝑦 represents the target variables including atomization energy, force, orbital energies. ℒ 𝐿 2 is the L2 regularization term on the network and ∆𝑷 is the deviation of SEQM parameters, a,b,c are the weights for these loss terms, we use a=1.0, b=10 -6 , and c=10 to have a stable training. 𝑓 and 𝑔 are the error functions, here we use root mean square error (RMSE) and mean absolute error (MAE) in f, and Mean Square Error (MSE) in g. As the molecule configuration is used as the input for HIPNN to generate SEQM parameters 𝑷 , there is an implicit additional term presented in the force calculation: 𝑭⃗⃗ = 𝜕 𝐸 𝑡𝑜𝑡 (𝑹⃗⃗ ;𝑷 ) 𝜕 𝑹⃗⃗ +( 𝜕 𝑷 𝜕 𝑹⃗⃗ ) 𝑇 𝜕 𝐸 𝑡𝑜𝑡 (𝑹⃗⃗ ;𝑷 ) 𝜕 𝑷 (7.3) The first term is the normal one for computing force, and the second one is coming from the backpropagation of the HIPNN network. With this full gradient of total energy on coordinates, the total energy is conserved during the dynamic simulations. 7.3 Training Details The training of the models closely follows the procedure detailed in ref 274 . For comparison of accuracy and extensibility, we train hybrid models (HIPNN+SEQM) and pure HIPNN model with different hyperparameters with a dataset generated with Density Functional Theory (DFT). The dataset consists of 618409 organic molecules with 5 to 18 atoms (1 to 13 carbon, nitrogen, oxygen atoms), as shown in Figure 7.2(a). The training dataset is generated with hybrid functional ωB97X with basis 6-31g* using Gaussian. 316 Due to the computational time cost, for the training of hybrid models, 10% molecules are randomly chosen (8% for training, 1% for validation and 1% for testing). As a tradeoff of GPU memory limit and learning efficiency, each batch has 256 molecules and it takes around 1 hour for each epoch on 4 NVIDIA GPU (CINT cluster on Darwin- Fe) for the hybrid models. The pure HIPNN models are trained with 10% and whole dataset for comparison. With the best trained models, to compare their performance, we apply them on 111 COMP6 datasets, which consists of much larger molecules as shown in Fig. 7.2 (b-d) and generate with the same DFT method. Figure 7.2: Distribution of molecule size used in training, testing. (a): The models are training on dataset with 5 to 18 HCNO atoms (Hydrogen, Carbon, Nitrogen, Oxygen), and 1-13 CHO atoms. (b-d): And the models are further tested on COMP6 datasets (ANI-MD, Drug Bank, GDB07to09, GDB10to13, s66x8, tripeptides), the 3 rows below show their distribution of total atoms (HCNO) and CNO atoms. (ANI-MD is not shown here as it has only 14 molecules, and GDB07to09, GDB10to13 are grouped together with s66x8 here in the second row as the molecule sizes are close) 7.4 Results and Discussion To examine the hybrid models, we first investigate the output dynamic SEQM parameters as shown in Figure 7.3. We apply the one of the hybrid models on COMP6 and extract the intermediate SEQM parameters and compare with the PM3 constant parameters. Other hybrid model has similar parameter distributions, for simplicity the parameter distributions from other models are not listed here. Among the 18 output SEQM parameters for each element (11 for Hydrogen), most of them are around the corresponding the PM3 constant parameters, as shown in Fig. 7.3(a-l). The deviations are within 30%, most of them are within 10% with respect to the PM3 constant ones. The distributions of these output parameters are asymmetric and not gaussian-like. 112 Most importantly, most of the distributions are multimodal, especially for parameters of Nitrogen as shown in Fig. 7.3(c, g, k), which shows the complex chemical environments in the compounds. Apart from this, some parameters like gss shown in Fig. 7.3(e-h) are systematically smaller than the corresponding PM3 constant ones, indicating the original PM3 parameters may have not been sufficiently optimized. Despite these relatively normal distributions, the exchange coulomb integral parameter hsp has a peak distribution around 0. During the training, hsp is converging to 0, and leads to instability in some numerical computation. As such, we put a hard lower bound 0.01 eV for hsp, this leads to the peak around 0.01 as shown in Fig. 7.3(m-o). hsp is used to construct several coulomb interaction integrals between s- and p-orbitals centered on different atoms, the vanishing of hsp shows that some integrals in the SEQM methods may be ignored or removed. This shows one possible direction to modify and improve the SEQM models. 113 Figure 7.3: Output SEQM parameter distributions from one of the best hybrid models. Here each column is for one element, for left to right, Hydrogen (green), Carbon (red), Nitrogen (blue), Oxygen (yellow). Eac row is for one type parameter, from top to bottom Uss: on-site energy term for s-orbital (unit eV), gss: direct coulomb integral between s-orbitals (unit eV), α: exponent term in nuclear interaction (unit 1/bohr radius), hsp: exchange coulomb integral between s- and p-orbitals (unit eV, only for CNO). The vertical black dash lines in (a-l) are to represent constant PM3 parameters for each of them. The numbers in (m-o) shows the constant PM3 hsp parameters. We apply these three models: original PM3, pure HIPNN, and hybrid HIPNN+SEQM, on COMP6 to compare their accuracies on predicting atomization energies and atomic forces. The 114 metric for them is listed in Table 7.1. During the training, the pure HIPNN model performs better than the hybrid one. However, when applying on COMP6, the hybrid model is better than the other two. Compared to the original PM3 model, the accuracy of the hybrid model is enhanced by around 50% on predicting energies and forces, as shown in Table 7.1. The COMB6 contains much larger molecules than ones used for training, especially in the subsets Tripeptides and Drug Bank. Both the original PM3 and the hybrid models show small difference on accuracy across the COMP6, while the pure HIPNN neural network model have larger errors on Tripeptides and Drug Bank. This indicates the quantum chemistry models are promising on addressing the extensibility issue presented in traditional machine learning models. One abnormal part is that the accuracy of hybrid model on COMP6 is even better than on test set from training. This may be because we apply a tighter threshold (10 -6 Eh) along with a combination of SCF algorithms on the SCF procedure when applying the model on COMB6. A set of SCF algorithms including constant mixing, adaptive mixing and Pulay convergers are used during the application phase to guarantee the success of SCF procedure, while during the training phase, the not converged samples are omitted. A tighter threshold will result in less noise on the single- particle density matrix D and lead to higher accuracy on the predicted quantities. This also shows that we may relieve this threshold to reduce the SCF steps during the training and make it less time-consuming. 115 Atomization Energy [eV] Force [eV/Å] RMSE MAE PM3 HIPNN HIPNN+ SEQM PM3 HIPNN HIPNN+ SEQM Test - - 0.12 0.08 0.20 0.15 - - 0.24 0.13 0.32 0.19 GDB+s66x8 0.28 0.19 0.18 0.11 0.13 0.09 0.63 0.40 0.33 0.21 0.24 0.17 Tripeptides 0.29 0.20 0.22 0.16 0.19 0.12 0.67 0.42 0.53 0.32 0.52 0.22 Drug Bank 0.29 0.22 0.30 0.20 0.18 0.13 0.53 0.35 0.52 0.32 0.28 0.19 Table 7.1: Accuracy on predicting atomization energies (unit: eV) and atomic forces (unit: eV/Å) from applying the original PM3 model, pure HIPNN model and the hybrid HIPNN+SEQM model on the test set during training and CDG+s66x8, Tripeptides, Drug Bank sets from COMP6. The extensibility of the models is further examined in Figure 7.4, which shows the predicting error of these models with respect to the size of chemical systems in COMP6. First of all, the error from the hybrid HIPNN+SEQM model is lowest over all the system size, for energies and atomic forces as shown in Figure 7.4(a) and (b). The hybrid model has errors which is around half of the errors from pure PM3 model and the hybrid one has close error with pure HIPNN model when applying to small systems with up to 20 atoms, but for larger systems, the pure HIPNN model performs as poorly as pure PM3. This indicates the hybrid model has better extensibility when applied to systems outside the training set, and efforts need to pay to tune the pure neural network for improving their extensibility. In Figure 7.4 (a) and (b), it shows the models have larger error on small systems with fewer than 10 atoms. It is due to the few data points in this range. There are only few small molecule systems in COMP6. In total, there are around 300 molecule configurations with fewer than 10 atoms, which constitutes less than 0.3 % of the whole COMP6 dataset. Fewer data points for 116 comparing and estimating the predicting errors lead to larger fluctuation in the results. Given the errors in the same magnitude, it is acceptable. One more thing is that the errors of these model scales differently with respect to the system size. In the system size range 5-90 given in the Figure 7.4 (a), all three models have error going up more or less linearly with system size for predicting atomization energies. As atomization energies are extensive, larger system with similar atom compounds will have larger atomization energy. And in the limit, energy should go up linearly with system size. This intrinsic linear scaling of atomization energy leads to its error from models follow the similar trend. However, the error on the forces behave differently. The errors from PM3 and hybrid models have more or less a constant error in the system size range 30-90. This is because COMP6 is generated with molecule coordinates fluctuate with given temperatures, as such the force fluctuate around within a range, and doesn’t show scaling with system sizes, so with its predicting error. And this suggests testing the transferability of predicting atomic forces, we may apply the models on systems with a range of temperatures or out of equilibrium to have a larger range on force distribution. Compared to the weak scaling of PM3 and hybrid model on predicting forces, the HIPNN model has errors going up with system sizes for this. This shows the methods built with quantum chemistry capture the essence of predicting forces, while pure neural network model here is relatively poor on this. 117 Figure 7.4: (a) MAE for predicted atomization energies for molecules in COMP6 over molecule size for three models: original PM3: black dashed line, pure HIPNN: red solid line, hybrid HIPNN+SEQM: blue solid line. (b) MAE for predicted atomic forces for these three models. 7.5 Conclusions In summary, we have incorporated artificial neural network with a simplified quantum chemistry model SEQM. SEQM is popular for its efficiency and thus can be used to study large systems with reasonable accuracy. The incorporation with neural network rejuvenizes this out-of- date method and improve its accuracy by 50%. At the same time, this hybrid model shows better extensibility compared to the pure neural network HIPNN when applied on systems outside the training set. This shows the promising of connecting with quantum chemistry with machine learning and adapts the domain knowledge in quantum chemistry to address the interpretability and extensibility problems in machine learning. Apart from directly hooking SEQM with HIPNN, we can also replace part of SEQM to fit scenarios where SEQM performs poorly. It has been well known that SEQM is problematic for description of non-valance interaction like van der Waals interaction, Hydrogen bonds. 315 A series of methods have been developed to address this, most of them modify the nuclei interaction forms. This unknown function form can be replaced with neural network, to automatically learn the 118 implicit function. 315 Given the flexibility of SEQM, we can also take its intermediate output for Δ- Machine Learning models. 312 Furthermore, SEQM has already been extended for excited state systems for which it is more computational demanding. 317 It is reasonable to adapt these extensions with machine learning technique and help address the weakness of SEQM models and at the same time provides methods for efficiently modeling systems which are out of reach with other methods. 119 Chapter 8: Closing Remarks Theoretical modeling in materials science and chemistry provides a unique perspective on studying nanoscale systems. In this thesis, we provide an overview of several methods used in modeling dynamics in nanoscale systems, and detail the effort we make to advance this field and expand their applications. We start with two applications to illustrate the modeling procedures and point out advantages and disadvantages of these modeling methods and the potential improvements we can make. Our approach begins with the method development for modeling of carrier-carrier scattering in low dimensional materials. And then we adapt modern machine learning techniques for help the model developments. As the most critical problem in computational material science is to have models with great efficiency and high accuracy, researchers expect to carry out modeling and calculations with feasible computation resource and with quantum-level accuracy to have trustful and useful results. The incorporation of machine learning in our method development paves the way for addressing this and exploring material science. The application works in Chapter 2 and 3 point out the effort we can make for the advancing of this field. In Chapter 2, we investigate the liquid phase exfoliation of layered material MoS2 with non-equilibrate molecular dynamics, show the atomic mechanism from the onset of exfoliation to the stress releasing in the systems and characterize essential properties of the system along the dynamics. Most importantly, we calibrate the force field for the description of solvent and crystal interfacial interaction. Apart from the results, the work shows the need of a convenient and easy way for the force field development in using classical molecular dynamics. In Chapter 3, we employ the non-adiabatic molecular dynamics to study and compare the carrier dynamics in low dimensional material carbon nanotube and graphene nanoribbon. The result shows the rigidity of 120 the whole system can have a great impact on carrier relaxation. This work gives us insight information on nanoscale systems, points out that the reduction of dimensions in nanomaterials may lead to fruitful phenomena with different physical origins, and shows the opportunities to introduce techniques for analyzing and extracting implicit information and dependency in conjunction with the modeling. Motivated by these application works, we move to the method development for studying nanoscale systems. In Chapter 4, we introduce the carrier-carrier interaction in the non-adiabatic molecular dynamics scheme. The strong carrier scattering in low dimensional materials coming from the spatial confinement dominates the carrier dynamics, while its theoretical study is limited with proper description and efficient and accurate modeling. The model we develop here successfully captures the multi-particle relaxation process in quantum dot CdSe and reveals detail carrier relaxation timescales and pathways. In Chapter 5, we adapt the unsupervised machine learning technique into non-adiabatic molecular dynamics and extract the dependency of carrier relaxation on the lattice deformation and vibration. The new scheme can provide quantitative estimation and comparison of factors affecting the carrier dynamics. And in Chapter 6 and 7, we provide a hybrid model with quantum chemistry methods and artificial neural network techniques. As claimed before, the major drawback of computational modeling of nanoscale systems is the high demand of computation resources. We can sacrifice accuracy with approximation and keep the essence of physics in developing models, and this is what people do from past several decades to now when trying to extend the model to study large and realistic systems. The emerging of machine learning enables us to push forward this boundary and have high accuracy and low cost in our models. The hybrid model provided in Chapter 6 and 7 gives us a demonstration example on organic systems and shows the power for modeling very large chemical systems with quantum 121 accuracy. And at the same time, further extension and improvement can be built on top of this for studying excited state phenomena like photovoltaics, chemical reactions, energy transfers for help material design and discovery. 225, 318 122 References 1. Allen, M. P.; Tildesley, D. J., Computer simulation of liquids. Oxford university press: 2017. 2. Zhou, G.; Rajak, P.; Susarla, S.; Ajayan, P. M.; Kalia, R. K.; Nakano, A.; Vashishta, P. Molecular Simulation of MoS2 Exfoliation. Sci. Rep. 2018, 8. 3. Hong, S.; Sheng, C.; Krishnamoorthy, A.; Rajak, P.; Tiwari, S.; Nomura, K.-i.; Misawa, M.; Shimojo, F.; Kalia, R. K.; Nakano, A. Chemical vapor deposition synthesis of MoS2 layers from the direct sulfidation of MoO3 surfaces using reactive molecular dynamics simulations. J. Phys. Chem. C 2018, 122, 7494-7503. 4. He, Y.; Nomura, K.-i.; Kalia, R. K.; Nakano, A.; Vashishta, P. Structure and dynamics of water confined in nanoporous carbon. Phys. Rev. Mater. 2018, 2, 115605. 5. Krishnamoorthy, A.; Rajak, P.; Norouzzadeh, P.; Singh, D. J.; Kalia, R. K.; Nakano, A.; Vashishta, P. Thermal conductivity of MoS2 monolayers from molecular dynamics simulations. AIP Advances 2019, 9, 035042. 6. Rajak, P.; Krishnamoorthy, A.; Nakano, A.; Vashishta, P.; Kalia, R. Structural phase transitions in a MoWSe2 monolayer: Molecular dynamics simulations and variational autoencoder analysis. Phys. Rev. B 2019, 100, 014108. 7. Marx, D.; Hutter, J., Ab initio molecular dynamics: basic theory and advanced methods. Cambridge University Press: 2009. 8. Martinez, T. J.; Ben-Nun, M.; Levine, R. Multi-electronic-state molecular dynamics: A wave function approach with applications. J. Phys. Chem. 1996, 100, 7884-7895. 9. White, A.; Tretiak, S.; Mozyrsky, D. Coupled wave-packets for non-adiabatic molecular dynamics: a generalization of Gaussian wave-packet dynamics to multiple potential energy surfaces. Chem. Sci. 2016, 7, 4905-4911. 10. Tully, J. C. Molecular dynamics with electronic transitions. J. Chem. Phys. 1990, 93, 1061- 1071. 11. Jaeger, H. M.; Fischer, S.; Prezhdo, O. V. Decoherence-induced surface hopping. J. Chem. Phys. 2012, 137, 22A545. 12. Nijjar, P.; Jankowska, J.; Prezhdo, O. V. Ehrenfest and classical path dynamics with decoherence and detailed balance. J. Chem. Phys. 2019, 150, 204124. 13. Wang, L.; Trivedi, D.; Prezhdo, O. V. Global flux surface hopping approach for mixed quantum-classical dynamics. J. Chem. Theory Comput. 2014, 10, 3598-3605. 14. Wang, L.; Beljonne, D. Flexible surface hopping approach to model the crossover from hopping to band-like transport in organic crystals. J. Phys. Chem. Lett. 2013, 4, 1888-1894. 15. Akimov, A. V.; Prezhdo, O. V. The PYXAID Program for Non-Adiabatic Molecular Dynamics in Condensed Matter Systems. J. Chem. Theory Comput. 2013, 9, 4959-4972. 16. Smith, J. S.; Isayev, O.; Roitberg, A. E. ANI-1: an extensible neural network potential with DFT accuracy at force field computational cost. Chem. Sci. 2017, 8, 3192-3203. 17. Sifain, A. E.; Lubbers, N.; Nebgen, B. T.; Smith, J. S.; Lokhov, A. Y.; Isayev, O.; Roitberg, A. E.; Barros, K.; Tretiak, S. Discovering a transferable charge assignment model using machine learning. J. Phys. Chem. Lett. 2018, 9, 4495-4501. 18. Zhou, G.; Chu, W.; Prezhdo, O. V. Structure Deformation Controls Charge Losses in MAPbI3: Unsupervised Machine Learning of Nonadiabatic Molecular Dynamics. ACS Ener. Lett. 2020. 123 19. Andersen, H. C. Rattle: A “velocity” version of the shake algorithm for molecular dynamics calculations. J. Comput. Phys. 1983, 52, 24-34. 20. Goldhirsch, I.; Sulem, P.-L.; Orszag, S. A. Stability and Lyapunov stability of dynamical systems: A differential approach and a numerical method. Physica D: Nonlinear Phenomena 1987, 27, 311-337. 21. Volz, S. G.; Chen, G. Molecular-dynamics simulation of thermal conductivity of silicon crystals. Phys. Rev. B 2000, 61, 2651. 22. Niu, H.; Piaggi, P. M.; Invernizzi, M.; Parrinello, M. Molecular dynamics simulations of liquid silica crystallization. Proc. Nat;. Acad. of Sci. 2018, 115, 5348-5352. 23. Neogi, A.; Mitra, N. Shock induced phase transition of water: Molecular dynamics investigation. Phys. Fluids 2016, 28, 027104. 24. Alavi, A.; Kohanoff, J.; Parrinello, M.; Frenkel, D. Ab initio molecular dynamics with excited electrons. Phys. Rev. Lett. 1994, 73, 2599. 25. Parandekar, P. V.; Tully, J. C. Mixed quantum-classical equilibrium. J. Chem. Phys. 2005, 122, 094102. 26. Wang, L. J.; Akimov, A.; Prezhdo, O. V. Recent Progress in Surface Hopping: 2011-2015. J. Phys. Chem. Lett. 2016, 7, 2100-2112. 27. Wu, S.; Kondo, Y.; Kakimoto, M.-a.; Yang, B.; Yamada, H.; Kuwajima, I.; Lambard, G.; Hongo, K.; Xu, Y.; Shiomi, J. Machine-learning-assisted discovery of polymers with high thermal conductivity using a molecular design algorithm. NPJ Comput. Mater. 2019, 5, 1-11. 28. Smith, J. S.; Nebgen, B. T.; Zubatyuk, R.; Lubbers, N.; Devereux, C.; Barros, K.; Tretiak, S.; Isayev, O.; Roitberg, A. E. Approaching coupled cluster accuracy with a general-purpose neural network potential through transfer learning. Nat. Commun. 2019, 10, 2903. 29. Nebgen, B.; Lubbers, N.; Smith, J. S.; Sifain, A. E.; Lokhov, A.; Isayev, O.; Roitberg, A. E.; Barros, K.; Tretiak, S. Transferable dynamic molecular charge assignment using deep neural networks. J. Chem. Theory Comput. 2018, 14, 4687-4698. 30. Engel, E.; Dreizler, R. M., Density Functional Theory. Springer: Berlin: 2013. 31. Senftle, T. P.; Hong, S.; Islam, M. M.; Kylasa, S. B.; Zheng, Y.; Shin, Y. K.; Junkermeier, C.; Engel-Herbert, R.; Janik, M. J.; Aktulga, H. M. The ReaxFF reactive force-field: development, applications and future directions. NPJ Comput. Mater. 2016, 2, 1-14. 32. Nicolosi, V.; Chhowalla, M.; Kanatzidis, M. G.; Strano, M. S.; Coleman, J. N. Liquid Exfoliation of Layered Materials. Science 2013, 340, 6139. 33. Han, J. T.; Jang, J. I.; Kim, H.; Hwang, J. Y.; Yoo, H. K.; Woo, J. S.; Choi, S.; Kim, H. Y.; Jeong, H. J.; Jeong, S. Y. Extremely efficient liquid exfoliation and dispersion of layered materials by unusual acoustic cavitation. Sci. Rep. 2014, 4, 5133. 34. Paton, K. R.; Varrla, E.; Backes, C.; Smith, R. J.; Khan, U.; O'Neill, A.; Boland, C.; Lotya, M.; Istrate, O. M.; King, P.; Higgins, T.; Barwich, S.; May, P.; Puczkarski, P.; Ahmed, I.; Moebius, M.; Pettersson, H.; Long, E.; Coelho, J.; O'Brien, S. E.; McGuire, E. K.; Sanchez, B. M.; Duesberg, G. S.; McEvoy, N.; Pennycook, T. J.; Downing, C.; Crossley, A.; Nicolosi, V.; Coleman, J. N. Scalable production of large quantities of defect-free few-layer graphene by shear exfoliation in liquids. Nat. Mater. 2014, 13, 624-630. 35. Niu, L. Y.; Coleman, J. N.; Zhang, H.; Shin, H.; Chhowalla, M.; Zheng, Z. J. Production of Two-Dimensional Nanomaterials via Liquid-Based Direct Exfoliation. Small 2016, 12, 272-293. 36. Ma, R.; Sasaki, T. Two-dimensional oxide and hydroxide nanosheets: controllable high- quality exfoliation, molecular assembly, and exploration of functionality. Acc. Chem. Res. 2015, 48, 136-143. 124 37. Coleman, J. N.; Lotya, M.; O’Neill, A.; Bergin, S. D.; King, P. J.; Khan, U.; Young, K.; Gaucher, A.; De, S.; Smith, R. J. Two-dimensional nanosheets produced by liquid exfoliation of layered materials. Science 2011, 331, 568-571. 38. Novoselov, K. S.; Fal, V.; Colombo, L.; Gellert, P.; Schwab, M.; Kim, K. A roadmap for graphene. Nature 2012, 490, 192-200. 39. Geim, A. K.; Novoselov, K. S., The rise of graphene. Nanoscience and Technology: a collection of reviews from nature journals, 2010, 11-19. 40. Huang, X.; Yin, Z.; Wu, S.; Qi, X.; He, Q.; Zhang, Q.; Yan, Q.; Boey, F.; Zhang, H. Graphene ‐based materials: synthesis, characterization, properties, and applications. Small 2011, 7, 1876-1902. 41. Britnell, L.; Gorbachev, R.; Jalil, R.; Belle, B.; Schedin, F.; Mishchenko, A.; Georgiou, T.; Katsnelson, M.; Eaves, L.; Morozov, S. Field-effect tunneling transistor based on vertical graphene heterostructures. Science 2012, 335, 947-950. 42. Georgiou, T.; Jalil, R.; Belle, B. D.; Britnell, L.; Gorbachev, R. V.; Morozov, S. V.; Kim, Y.-J.; Gholinia, A.; Haigh, S. J.; Makarovsky, O. Vertical field-effect transistor based on graphene–WS 2 heterostructures for flexible and transparent electronics. Nat. Nanotechnol. 2013, 8, 100. 43. Kafy, A.; Sadasivuni, K. K.; Kim, H.-C.; Akther, A.; Kim, J. Designing flexible energy and memory storage materials using cellulose modified graphene oxide nanocomposites. Phys. Chem. Chem. Phys. 2015, 17, 5923-5931. 44. Lee, K. H.; Shin, H.-J.; Lee, J.; Lee, I.-y.; Kim, G.-H.; Choi, J.-Y.; Kim, S.-W. Large-scale synthesis of high-quality hexagonal boron nitride nanosheets for large-area graphene electronics. Nano Lett. 2012, 12, 714-718. 45. Le, L. T.; Ervin, M. H.; Qiu, H.; Fuchs, B. E.; Lee, W. Y. Graphene supercapacitor electrodes fabricated by inkjet printing and thermal reduction of graphene oxide. Electrochem. Commun. 2011, 13, 355-358. 46. Shen, J.; He, Y.; Wu, J.; Gao, C.; Keyshar, K.; Zhang, X.; Yang, Y.; Ye, M.; Vajtai, R.; Lou, J. Liquid phase exfoliation of two-dimensional materials by directly probing and matching surface tension components. Nano Lett. 2015, 15, 5449-5454. 47. Wang, M.; Xu, X.; Ge, Y.; Dong, P.; Baines, R.; Ajayan, P. M.; Ye, M.; Shen, J. Surface tension components ratio: an efficient parameter for direct liquid phase exfoliation. ACS Appl. Mater. Interfaces 2017, 9, 9168-9175. 48. Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 1995, 117, 1-19. 49. Plimpton, S.; Thompson, A.; Crozier, P.; Kohlmeyer, A. LAMMPS molecular dynamics simulator. URL http://lammps. sandia. gov 2011. 50. Rybakov, A.; Rybakov, I. Polymorphism of shocked water. Eur. J. Mech. B, Fluids 1995, 14, 323-332. 51. Gleason, A.; Bolme, C. A.; Galtier, E.; Lee, H. J.; Granados, E.; Dolan, D.; Seagle, C.; Ao, T.; Ali, S.; Lazicki, A. Compression freezing kinetics of water to ice VII. Phys. Rev. Lett. 2017, 119, 025701. 52. Vedadi, M.; Choubey, A.; Nomura, K.-I.; Kalia, R.; Nakano, A.; Vashishta, P.; Van Duin, A. Structure and dynamics of shock-induced nanobubble collapse in water. Phys. Rev. Lett. 2010, 105, 014503. 53. Kodama, T.; Tomita, Y. Cavitation bubble behavior and bubble–shock wave interaction near a gelatin surface as a study of in vivo bubble dynamics. Appl. Phys. B 2000, 70, 139-149. 125 54. Ohl, C.; Ikink, R. Shock-wave-induced jetting of micron-size bubbles. Phys. Rev. Lett. 2003, 90, 214502. 55. Ohl, C.-D.; Arora, M.; Dijkink, R.; Janve, V.; Lohse, D. Surface cleaning from laser- induced cavitation bubbles. Appl. Phys. Lett. 2006, 89, 074102. 56. Suslick, K. S.; Flannigan, D. J. Inside a collapsing bubble: sonoluminescence and the conditions during cavitation. Annu. Rev. Phys. Chem. 2008, 59, 659-683. 57. Abascal, J. L.; Vega, C. A general purpose model for the condensed phases of water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. 58. Liang, T.; Phillpot, S. R.; Sinnott, S. B. Parametrization of a reactive many-body potential for Mo–S systems. Phys. Rev. B 2009, 79, 245110. 59. Stewart, J. A.; Spearot, D. Atomistic simulations of nanoindentation on the basal plane of crystalline molybdenum disulfide (MoS2). Modell. Simul. Mater. Sci. Eng. 2013, 21, 045003. 60. Jorgensen, W., OPLS All-Atom Parameters for Organic Molecules, Ions, Peptides & Nucleic Acids. Yale University, New Haven, CT: 2009. 61. Luan, B.; Zhou, R. Wettability and friction of water on a MoS2 nanosheet. Appl. Phys. Lett. 2016, 108, 131601. 62. Yeh, I.-C.; Berkowitz, M. L. Ewald summation for systems with slab geometry. J. Chem. Phys. 1999, 111, 3155-3162. 63. Halim, U.; Zheng, C. R.; Chen, Y.; Lin, Z.; Jiang, S.; Cheng, R.; Huang, Y.; Duan, X. A rational design of cosolvent exfoliation of layered materials by directly probing liquid–solid interaction. Nat. Commun. 2013, 4, 1-7. 64. Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO–the Open Visualization Tool. Modell. Simul. Mater. Sci. Eng. 2009, 18, 015012. 65. Ahrens, J.; Geveci, B.; Law, C. Paraview: An end-user tool for large data visualization. The visualization handbook 2005, 717. 66. Stukowski, A. Computational analysis methods in atomistic modeling of crystals. JOM 2014, 66, 399-407. 67. Barber, C. B.; Dobkin, D. P.; Huhdanpaa, H. The quickhull algorithm for convex hulls. ACM Trans. Math. Software (TOMS) 1996, 22, 469-483. 68. Shimizu, F.; Ogata, S.; Li, J. Theory of shear banding in metallic glasses and molecular dynamics calculations. Mater. Trans. 2007, 0710160231-0710160231. 69. Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; Jiang, D.; Zhang, Y.; Dubonos, S. V.; Grigorieva, I. V.; Firsov, A. A. Electric field effect in atomically thin carbon films. Science 2004, 306, 666-669. 70. Yu, D. S.; Dai, L. M. Self-Assembled Graphene/Carbon Nanotube Hybrid Films for Supercapacitors. J. Phys. Chem. Lett. 2010, 1, 467-470. 71. Castro Neto, A. H.; Guinea, F.; Peres, N. M. R.; Novoselov, K. S.; Geim, A. K. The electronic properties of graphene. Rev. Mod. Phys. 2009, 81, 109-162. 72. Xiang, Q. J.; Yu, J. G. Graphene-Based Photocatalysts for Hydrogen Generation. J. Phys. Chem. Lett. 2013, 4, 753-759. 73. Iijima, S.; Ichihashi, T. Single-shell Carbon Nanotubes of 1-nm Diameter. Nature 1993, 363, 603-605. 74. Bolotin, K. I.; Sikes, K. J.; Jiang, Z.; Klima, M.; Fudenberg, G.; Hone, J.; Kim, P.; Stormer, H. Ultrahigh electron mobility in suspended graphene. Solid State Commun. 2008, 146, 351-355. 126 75. Bachilo, S. M.; Strano, M. S.; Kittrell, C.; Hauge, R. H.; Smalley, R. E.; Weisman, R. B. Structure-assigned optical spectra of single-walled carbon nanotubes. Science 2002, 298, 2361- 2366. 76. Novoselov, K. S.; Geim, A. K.; Morozov, S.; Jiang, D.; Katsnelson, M. I.; Grigorieva, I.; Dubonos, S.; Firsov; AA Two-dimensional gas of massless Dirac fermions in graphene. Nature 2005, 438, 197. 77. Kane, C. L.; Mele, E. J. Quantum spin Hall effect in graphene. Phys. Rev. Lett. 2005, 95, 226801. 78. Dai, H. J. Carbon nanotubes: Synthesis, integration, and properties. Acc. Chem. Res. 2002, 35, 1035-1044. 79. Zhang, Y.; Tan, Y.-W.; Stormer, H. L.; Kim, P. Experimental observation of the quantum Hall effect and Berry's phase in graphene. Nature 2005, 438, 201. 80. McEuen, P. L.; Fuhrer, M. S.; Park, H. Single-walled carbon nanotube electronics. IEEE Trans. Nanotechnol. 2002, 1, 78-85. 81. Park, S.; Vosguerichian, M.; Bao, Z. A review of fabrication and applications of carbon nanotube film-based flexible electronics. Nanoscale 2013, 5, 1727-1752. 82. Rakhi, R., Preparation and properties of manipulated carbon nanotube composites and applications. Nanocarbon Its Compos., Elsevier: 2019, 489-520. 83. Allen, M. J.; Tung, V. C.; Kaner, R. B. Honeycomb carbon: a review of graphene. Chem. Rev. 2009, 110, 132-145. 84. Jia, X.; Campos-Delgado, J.; Terrones, M.; Meunier, V.; Dresselhaus, M. S. Graphene edges: a review of their fabrication and characterization. Nanoscale 2011, 3, 86-95. 85. Chen, Y.-C.; De Oteyza, D. G.; Pedramrazi, Z.; Chen, C.; Fischer, F. R.; Crommie, M. F. Tuning the band gap of graphene nanoribbons synthesized from molecular precursors. ACS Nano 2013, 7, 6123-6128. 86. Strano, M. S.; Dyke, C. A.; Usrey, M. L.; Barone, P. W.; Allen, M. J.; Shan, H.; Kittrell, C.; Hauge, R. H.; Tour, J. M.; Smalley, R. E. Electronic structure control of single-walled carbon nanotube functionalization. Science 2003, 301, 1519-1522. 87. Sun, D. M.; Liu, C.; Ren, W. C.; Cheng, H. M. A review of carbon nanotube ‐and graphene ‐based flexible thin ‐film transistors. Small 2013, 9, 1188-1205. 88. Lee, S.-H.; Sridhar, V.; Jung, J.-H.; Karthikeyan, K.; Lee, Y.-S.; Mukherjee, R.; Koratkar, N.; Oh, I.-K. Graphene–nanotube–iron hierarchical nanostructure as lithium ion battery anode. ACS Nano 2013, 7, 4242-4251. 89. Zhao, M.-Q.; Liu, X.-F.; Zhang, Q.; Tian, G.-L.; Huang, J.-Q.; Zhu, W.; Wei, F. Graphene/single-walled carbon nanotube hybrids: one-step catalytic growth and applications for high-rate Li–S batteries. ACS Nano 2012, 6, 10759-10769. 90. Lin, J.; Peng, Z.; Xiang, C.; Ruan, G.; Yan, Z.; Natelson, D.; Tour, J. M. Graphene nanoribbon and nanostructured SnO2 composite anodes for lithium ion batteries. ACS Nano 2013, 7, 6001-6006. 91. Dillon, A. C.; Jones, K.; Bekkedahl, T.; Kiang, C.; Bethune, D.; Heben, M. Storage of hydrogen in single-walled carbon nanotubes. Nature 1997, 386, 377. 92. Lee, H.; Ihm, J.; Cohen, M. L.; Louie, S. G. Calcium-decorated graphene-based nanostructures for hydrogen storage. Nano Lett. 2010, 10, 793-798. 93. Jensen, S. A.; Ulbricht, R.; Narita, A.; Feng, X.; Mü llen, K.; Hertel, T.; Turchinovich, D.; Bonn, M. Ultrafast photoconductivity of graphene nanoribbons and carbon nanotubes. Nano Lett. 2013, 13, 5925-5930. 127 94. Zhu, Z.; Crochet, J.; Arnold, M. S.; Hersam, M. C.; Ulbricht, H.; Resasco, D.; Hertel, T. Pump-probe spectroscopy of exciton dynamics in (6, 5) carbon nanotubes. J. Phys. Chem. C 2007, 111, 3831-3835. 95. Hertel, T.; Himmelein, S.; Ackermann, T.; Stich, D.; Crochet, J. Diffusion Limited Photoluminescence Quantum Yields in 1-D Semiconductors: Single-Wall Carbon Nanotubes. ACS Nano 2010, 4, 7161-7168. 96. Kim, Y.; Velizhanin, K. A.; He, X. W.; Sarpkaya, I.; Yomogida, Y.; Tanaka, T.; Kataura, H.; Doorn, S. K.; Htoon, H. Photoluminescence Intensity Fluctuations and Temperature- Dependent Decay Dynamics of Individual Carbon Nanotube sp3 Defects. J. Phys. Chem. Lett. 2019, 10, 1423-1430. 97. del Valle, D. G. F.; Moretti, L.; Maqueira-Albo, I.; Aluicio-Sarduy, E.; Kriegel, I.; Lanzani, G.; Scotognella, F. Ultrafast Hole Transfer from (6,5) SWCNT to P3HT:PCBM Blend by Resonant Excitation. J. Phys. Chem. Lett. 2016, 7, 3353-3358. 98. Mann, C.; Hertel, T. 13 nm Exciton Size in (6,5) Single-Wall Carbon Nanotubes. J. Phys. Chem. Lett. 2016, 7, 2276-2280. 99. Castillo, M.; Pho, C.; Naumov, A. V.; Dzyuba, S. V. Modulating Chirality-Selective Photoluminescence of Single-Walled Carbon Nanotubes by Ionic Liquids. J. Phys. Chem. Lett. 2018, 9, 6689-6694. 100. Kryjevski, A.; Mihaylov, D.; Kilin, D. Dynamics of Charge Transfer and Multiple Exciton Generation in the Doped Silicon Quantum Dot-Carbon Nanotube System: Density Functional Theory-Based Computation. J. Phys. Chem. Lett. 2018, 9, 5759-5764. 101. Habenicht, B. F.; Prezhdo, O. V. Nonradiative quenching of fluorescence in a semiconducting carbon nanotube: A time-domain ab initio study. Phys. Rev. Lett. 2008, 100. 102. Chaban, V. V.; Pal, S.; Prezhdo, O. V. Laser-Induced Explosion of Nitrated Carbon Nanotubes: Nonadiabatic and Reactive Molecular Dynamics Simulations. J. Am. Chem. Soc. 2016, 138, 15927-15934. 103. Lauret, J. S.; Voisin, C.; Cassabois, G.; Delalande, C.; Roussignol, P.; Jost, O.; Capes, L. Ultrafast carrier dynamics in single-wall carbon nanotubes. Phys. Rev. Lett. 2003, 90. 104. Pal, S.; Casanova, D.; Prezhdo, O. V. Effect of Aspect Ratio on Multiparticle Auger Recombination in Single-Walled Carbon Nanotubes: Time Domain Atomistic Simulation. Nano Lett. 2018, 18, 58-63. 105. Ulbricht, R.; Hendry, E.; Shan, J.; Heinz, T. F.; Bonn, M. Carrier dynamics in semiconductors studied with time-resolved terahertz spectroscopy. Rev. Mod. Phys. 2011, 83, 543- 586. 106. Pal, S.; Trivedi, D. J.; Akimov, A. V.; Aradi, B.; Frauenheim, T.; Prezhdo, O. V. Nonadiabatic Molecular Dynamics for Thousand Atom Systems: A Tight-Binding Approach toward PYXAID. J. Chem. Theory Comput. 2016, 12, 1436-1448. 107. Sarkar, R.; Habib, M.; Pal, S.; Prezhdo, O. V. Ultrafast, asymmetric charge transfer and slow charge recombination in porphyrin/CNT composites demonstrated by time-domain atomistic simulation. Nanoscale 2018, 10, 12683-12694. 108. Xiong, W. J.; Du, L. L.; Lo, K. C.; Shi, H. T.; Takaya, T.; Iwata, K.; Chan, W. K.; Phillips, D. L. Control of Electron Flow Direction in Photoexcited Cycloplatinated Complex Containing Conjugated Polymer-Single-Walled Carbon Nanotube Hybrids. J. Phys. Chem. Lett. 2018, 9, 3819-3824. 128 109. Long, R.; Casanova, D.; Fang, W. H.; Prezhdo, O. V. Donor Acceptor Interaction Determines the Mechanism of Photoinduced Electron Injection from Graphene Quantum Dots into TiO2:pi-Stacking Supersedes Covalent Bonding. J. Am. Chem. Soc. 2017, 139, 2619-2629. 110. Habenicht, B. F.; Craig, C. F.; Prezhdo, O. V. Time-domain ab initio simulation of electron and hole relaxation dynamics in a single-wall semiconducting carbon nanotube. Phys. Rev. Lett. 2006, 96, 187401. 111. Hertel, T.; Moos, G. Electron-phonon interaction in single-wall carbon nanotubes: A time- domain study. Phys. Rev. Lett. 2000, 84, 5002. 112. Tully, J. C.; Preston, R. K. Trajectory Surface Hopping Approach to Nonadiabatic Molecular Collisions: Reaction of H + with D2. J. Chem. Phys. 1971, 55, 562-&. 113. Bittner, E. R.; Rossky, P. J. Quantum decoherence in mixed quantum ‐classical systems: Nonadiabatic processes. J. Chem. Phys. 1995, 103, 8130-8143. 114. Prezhdo, O. V.; Rossky, P. J. Mean-field molecular dynamics with surface hopping. J. Chem. Phys. 1997, 107, 825-834. 115. Prezhdo, O. V. Mean field approximation for the stochastic Schrodinger equation. J. Chem. Phys. 1999, 111, 8366-8377. 116. Akimov, A. V.; Long, R.; Prezhdo, O. V. Coherence penalty functional: A simple method for adding decoherence in Ehrenfest dynamics. J. Chem. Phys. 2014, 140. 117. Brooksby, C.; Prezhdo, O. V. Quantized mean-field approximation. Chem. Phys. Lett. 2001, 346, 463-469. 118. Craig, C. F.; Duncan, W. R.; Prezhdo, O. V. Trajectory surface hopping in the time- dependent Kohn-Sham approach for electron-nuclear dynamics. Phys. Rev. Lett. 2005, 95. 119. Akimov, A. V.; Prezhdo, O. V. Advanced Capabilities of the PYXAID Program: Integration Schemes, Decoherenc:e Effects, Multiexcitonic States, and Field-Matter Interaction. J. Chem. Theory Comput. 2014, 10, 789-804. 120. Kresse, G.; Furthmuller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B 1996, 54, 11169-11186. 121. Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865-3868. 122. Blochl, P. E. Projector augmented-wave method. Phys. Rev. B 1994, 50, 17953-17979. 123. Akimov, A. V.; Prezhdo, O. V. Persistent Electronic Coherence Despite Rapid Loss of Electron-Nuclear Correlation. J. Phys. Chem. Lett. 2013, 4, 3857-3864. 124. Li, L. Q.; Long, R.; Prezhdo, O. V. Why Chemical Vapor Deposition Grown MoS2 Samples Outperform Physical Vapor Deposition Samples: Time-Domain ab Initio Analysis. Nano Lett. 2018, 18, 4008-4014. 125. He, J. L.; Vasenko, A. S.; Long, R.; Prezhdo, O. V. Halide Composition Controls Electron- Hole Recombination in Cesium-Lead Halide Perovskite Quantum Dots: A Time Domain Ab lnitio Study. J. Phys. Chem. Lett. 2018, 9, 1872-1879. 126. Zhang, L. L.; Vasenko, A. S.; Zhao, J.; Prezhdo, O. V. Mono-Elemental Properties of 2D Black Phosphorus Ensure Extended Charge Carrier Lifetimes under Oxidation: Time-Domain Ab Initio Analysis. J. Phys. Chem. Lett. 2019, 10, 1083-1091. 127. Jankowska, J.; Prezhdo, O. V. Real-Time Atomistic Dynamics of Energy Flow in an STM Setup: Revealing the Mechanism of Current-Induced Molecular Emission. J. Phys. Chem. Lett. 2018, 9, 3591-3597. 129 128. Li, W.; Long, R.; Hou, Z. F.; Tang, J. F.; Prezhdo, O. V. Influence of Encapsulated Water on Luminescence Energy, Line Width, and Lifetime of Carbon Nanotubes: Time Domain Ab Initio Analysis. J. Phys. Chem. Lett. 2018, 9, 4006-4013. 129. Li, W.; Long, R.; Tang, J.; Prezhdo, O. V. Influence of Defects on Excited-State Dynamics in Lead Halide Perovskites: Time-Domain ab Initio Studies. J. Phys. Chem. Lett. 2019, 10, 3788- 3804. 130. Wang, Y. T.; Fang, W. H.; Long, R.; Prezhdo, O. V. Symmetry Breaking at MAPbI(3) Perovskite Grain Boundaries Suppresses Charge Recombination: Time-Domain ab Initio Analysis. J. Phys. Chem. Lett. 2019, 10, 1617-1623. 131. Liu, L. H.; Fang, W. H.; Long, R.; Prezhdo, O. V. Lewis Base Passivation of Hybrid Halide Perovskites Slows Electron-Hole Recombination: Time-Domain Ab lnitio Analysis. J. Phys. Chem. Lett. 2018, 9, 1164-1171. 132. Zhang, Z. S.; Liu, L. H.; Fang, W. H.; Long, R.; Tokina, M. V.; Prezhdo, O. V. Plasmon- Mediated Electron Injection from Au Nanorods into MoS2: Traditional versus Photoexcitation Mechanism. Chem 2018, 4, 1112-1127. 133. Fischer, S. A.; Duncan, W. R.; Prezhdo, O. V. Ab Initio Nonadiabatic Molecular Dynamics of Wet-Electrons on the TiO2 Surface. J. Am. Chem. Soc. 2009, 131, 15483-15491. 134. Kilin, D. S.; Tsemekhman, K.; Prezhdo, O. V.; Zenkevich, E. I.; von Borczyskowski, C. Ab initio study of exciton transfer dynamics from a core-shell semiconductor quantum dot to a porphyrin-sensitizer. J. Photochem. Photobiol., A. 2007, 190, 342-351. 135. Wang, L. J.; Long, R.; Prezhdo, O. V., Time-Domain Ab Initio Modeling of Photoinduced Dynamics at Nanoscale Interfaces. Annu. Rev. Phys. Chem., 2015, 66, 549-579. 136. Akimov, A. V.; Asahi, R.; Jinnouchi, R.; Prezhdo, O. V. What Makes the Photocatalytic CO2 Reduction on N-Doped Ta2O5 Efficient: Insights from Nonadiabatic Molecular Dynamics. J. Am. Chem. Soc. 2015, 137, 11517-11525. 137. Kilina, S.; Kilin, D.; Tretiak, S. Light-Driven and Phonon-Assisted Dynamics in Organic and Semiconductor Nanostructures. Chem. Rev. 2015, 115, 5929-5978. 138. Nayyar, I. H.; Batista, E. R.; Tretiak, S.; Saxena, A.; Smith, D. L.; Martin, R. L. Localization of Electronic Excitations in Conjugated Polymers Studied by DFT. J. Phys. Chem. Lett. 2011, 2, 566-571. 139. Kilina, S. V.; Neukirch, A. J.; Habenicht, B. F.; Kilin, D. S.; Prezhdo, O. V. Quantum Zeno Effect Rationalizes the Phonon Bottleneck in Semiconductor Quantum Dots. Phys. Rev. Lett. 2013, 110. 140. Prezhdo, O. V. Quantum anti-zeno acceleration of a chemical reaction. Phys. Rev. Lett. 2000, 85, 4413. 141. Xue, Y. Q.; Datta, S.; Ratner, M. A. First-principles based matrix Green's function approach to molecular electronic devices: general formalism. Chem. Phys. 2002, 281, 151-170. 142. Efros, A. L.; Kharchenko, V. A.; Rosen, M. Breaking the phonon bottleneck in nanometer quantum dots: Role of Auger-like processes. Solid State Comm. 1995, 93, 281-284. 143. Klimov, V. I.; Mikhailovsky, A. A.; McBranch, D. W.; Leatherdale, C. A.; Bawendi, M. G. Quantization of multiparticle Auger rates in semiconductor quantum dots. Science 2000, 287, 1011-1013. 144. Cragg, G. E.; Efros, A. L. Suppression of Auger processes in confined structures. Nano Lett. 2010, 10, 313-317. 145. Philbin, J. P.; Rabani, E. Electron-Hole Correlations Govern Auger Recombination in Nanostructures. Nano Lett. 2018, 18, 7889-7895. 130 146. Covito, F.; Perfetto, E.; Rubio, A.; Stefanucci, G. Real-time dynamics of Auger wave packets and decays in ultrafast charge migration processes. Phys. Rev. A 2018, 97. 147. Gao, J. B.; Kidon, L.; Rabani, E.; Alivisatos, A. P. Ultrahigh Hot Carrier Transient Photocurrent in Nanocrystal Arrays by Auger Recombination. Nano Lett. 2019, 19, 4804-4810. 148. Nozik, A. J.; Beard, M. C.; Luther, J. M.; Law, M.; Ellingson, R. J.; Johnson, J. C. Semiconductor Quantum Dots and Quantum Dot Arrays and Applications of Multiple Exciton Generation to Third-Generation Photovoltaic Solar Cells. Chem. Rev. 2010, 110, 6873-6890. 149. Pandey, A.; Guyot-Sionnest, P. Slow Electron Cooling in Colloidal Quantum Dots. Science 2008, 322, 929-932. 150. Wang, L.; Chen, Z.; Liang, G.; Li, Y.; Lai, R.; Ding, T.; Wu, K. Observation of a phonon bottleneck in copper-doped colloidal quantum dots. Nat. Commun. 2019, 10, 1-8. 151. Trivedi, D. J.; Wang, L.; Prezhdo, O. V. Auger-mediated electron relaxation is robust to deep hole traps: time-domain ab initio study of CdSe quantum dots. Nano Lett. 2015, 15, 2086- 2091. 152. Mocker, M.; Lemke, F. Theoretical results on the efficiency of Auger recombination in low dimensional PbSnTe structures. Superlattices Microstruct. 1991, 10, 231-234. 153. Landau, L. D.; Lifshitz, E. M., Quantum mechanics: non-relativistic theory. Elsevier: 2013, 3. 154. Danilov, L.; Zegrya, G. Theoretical study of auger recombination processes in deep quantum wells. Semiconductors 2008, 42, 550-556. 155. Philbin, J. P.; Rabani, E. Auger Recombination Lifetime Scaling for Type I and Quasi- Type II Core/Shell Quantum Dots. J. Phys. Chem. Lett. 2020, 11, 5132-5138. 156. Hyeon-Deuk, K.; Prezhdo, O. V. Time-domain ab initio study of auger and phonon-assisted auger processes in a semiconductor quantum dot. Nano Lett. 2011, 11, 1845-1850. 157. Ben-Shahar, Y.; Philbin, J. P.; Scotognella, F.; Ganzer, L.; Cerullo, G.; Rabani, E.; Banin, U. Charge Carrier Dynamics in Photocatalytic Hybrid Semiconductor-Metal Nanorods: Crossover from Auger Recombination to Charge Transfer. Nano Lett. 2018, 18, 5211-5216. 158. Balan, A. D.; Eshet, H.; Olshansky, J. H.; Lee, Y. V.; Rabani, E.; Alivisatos, A. P. Effect of Thermal Fluctuations on the Radiative Rate in Core/Shell Quantum Dots. Nano Lett. 2017, 17, 1629-1636. 159. Zhang, Z. S.; Long, R.; Tokina, M. V.; Prezhdo, O. V. Interplay between Localized and Free Charge Carriers Can Explain Hot Fluorescence in the CH3NH3PbBr3 Perovskite: Time- Domain Ab Initio Analysis. J. Am. Chem. Soc. 2017, 139, 17327-17333. 160. Kilina, S.; Batista, E. R.; Yang, P.; Tretiak, S.; Saxena, A.; Martin, R. L.; Smith, D. L. Electronic structure of self-assembled amorphous polyfluorenes. ACS Nano 2008, 2, 1381-1388. 161. Hyeon-Deuk, K.; Prezhdo, O. V. Multiple exciton generation and recombination dynamics in small si and cdse quantum dots: An ab initio time-domain study. ACS Nano 2012, 6, 1239-1250. 162. Zhu, H.; Yang, Y.; Hyeon-Deuk, K.; Califano, M.; Song, N.; Wang, Y.; Zhang, W.; Prezhdo, O. V.; Lian, T. Auger-assisted electron transfer from photoexcited semiconductor quantum dots. Nano Lett. 2014, 14, 1263-1269. 163. Dong, S.; Pal, S.; Lian, J.; Chan, Y.; Prezhdo, O. V.; Loh, Z.-H. Sub-picosecond Auger- mediated hole-trapping dynamics in colloidal CdSe/CdS core/shell nanoplatelets. ACS Nano 2016, 10, 9370-9378. 164. Helbig, N.; Fuks, J. I.; Tokatly, I. V.; Appel, H.; Gross, E. K. U.; Rubio, A. Time-dependent density-functional and reduced density-matrix methods for few electrons: Exact versus adiabatic approximations. Chem. Phys. 2011, 391, 1-10. 131 165. Cooney, R. R.; Sewall, S. L.; Anderson, K. E.; Dias, E. A.; Kambhampati, P. Breaking the phonon bottleneck for holes in semiconductor quantum dots. Phys. Rev. Lett. 2007, 98, 177403. 166. Schaller, R. D.; Pietryga, J. M.; Goupalov, S. V.; Petruska, M. A.; Ivanov, S. A.; Klimov, V. I. Breaking the phonon bottleneck in semiconductor nanocrystals via multiphonon emission induced by intrinsic nonadiabatic interactions. Phys. Rev. Lett. 2005, 95, 196401. 167. Kohn, W.; Sham, L. J. Self-consistent equations including exchange and correlation effects. Phys. Rev. 1965, 140, A1133. 168. Kresse, G. Ab initio molecular dynamics for liquid metals. J. Non-Cryst. Solids 1995, 193, 222-229. 169. Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter. 2009, 21, 395502. 170. Wing, D.; Haber, J. B.; Noff, R.; Barker, B.; Egger, D. A.; Ramasubramaniam, A.; Louie, S. G.; Neaton, J. B.; Kronik, L. Comparing time-dependent density functional theory with many- body perturbation theory for semiconductors: Screened range-separated hybrids and the GW plus Bethe-Salpeter approach. Phys. Rev. Mater. 2019, 3. 171. Krylov, A. I. Equation-of-motion coupled-cluster methods for open-shell and electronically excited species: The Hitchhiker's guide to Fock space. Annu. Rev. Phys. Chem. 2008, 59, 433-462. 172. 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. 173. Castro, A.; Marques, M. A. L.; Alonso, J. A.; Rubio, A. Optical Properties of Nanostructures from Time-Dependent Density Functional Theory. J. Comput. Theor. Nanosci. 2004, 1, 231-255. 174. Marques, M. A. L.; Gross, E. K. U. Time-dependent density functional theory. Annu. Rev. Phys. Chem. 2004, 55, 427-455. 175. Marques, M. A. L.; Castro, A.; Bertsch, G. F.; Rubio, A. Octopus: a first-principles tool for excited electron-ion dynamics. Comput. Phys. Commun. 2003, 151, 60-78. 176. Tancogne-Dejean, N.; Oliveira, M. J. T.; Andrade, X.; Appel, H.; Borca, C. H.; Le Breton, G.; Buchholz, F.; Castro, A.; Corni, S.; Correa, A. A.; De Giovannini, U.; Delgado, A.; Bich, F. G.; Flick, J.; Gil, G.; Gomez, A.; Helbig, N.; Hubener, H.; Jestadt, R.; Jornet-Somoza, J.; Larsen, A. H.; Lebedeva, I. V.; Luders, M.; Marques, M. A. L.; Ohlmann, S. T.; Pipolo, S.; Rampp, M.; Rozzi, C. A.; Strubbe, D. A.; Sato, S. A.; Schafer, C.; Theophilou, I.; Welden, A.; Rubio, A. Octopus, a computational framework for exploring light-driven phenomena and quantum dynamics in extended and finite systems. J. Chem. Phys. 2020, 152. 177. Wang, L. J.; Prezhdo, O. V. A Simple Solution to the Trivial Crossing Problem in Surface Hopping. J. Phys. Chem. Lett. 2014, 5, 713-719. 178. Fernandez-Alberti, S.; Roitberg, A. E.; Nelson, T.; Tretiak, S. Identification of unavoided crossings in nonadiabatic photoexcited dynamics involving multiple electronic states in polyatomic conjugated molecules. J. Chem. Phys. 2012, 137. 179. Li, X. S.; Tully, J. C.; Schlegel, H. B.; Frisch, M. J. Ab initio Ehrenfest dynamics. J. Chem. Phys. 2005, 123. 180. Gisin, N.; Percival, I. C. The quantum-state diffusion model applied to open systems. J. Phys. A-Math. Gen. 1992, 25, 5677-5691. 132 181. Akimoy, A. V. A Simple Phase Correction Makes a Big Difference in Nonadiabatic Molecular Dynamics. J. Phys. Chem. Lett. 2018, 9, 6096-6102. 182. Kresse, G.; Furthmuller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. 1996, 6, 15-50. 183. Kresse, G.; Hafner, J. Ab initio molecular-dynamics simulation of the liquid-metal– amorphous-semiconductor transition in germanium. Phys. Rev. B 1994, 49, 14251. 184. Chu, W. B., Zheng Q., Akimov A. V., Zhao, J., Saidi, W., Prezhdo, O. V. Accurate Computation of Nonadiabatic Coupling with Projector Augmented-Wave Pseudopotentials. J. Phys. Chem. Lett. 2020, in press 185. Bystrom, K.; Broberg, D.; Dwaraknath, S.; Persson, K. A.; Asta, M. Pawpyseed: Perturbation-extrapolation band shifting corrections for point defect calculations. arXiv preprint arXiv:1904.11572 2019. 186. Zhang, X.; Lu, G. First-order nonadiabatic couplings in extended systems by time- dependent density functional theory. J. Chem. Phys. 2018, 149. 187. Isborn, C. M.; Prezhdo, O. V. Charging Quenches Multiple Exciton Generation in Semiconductor Nanocrystals: First-Principles Calculations on Small PbSe Clusterse. J. Phys. Chem. C 2009, 113, 12617-12621. 188. Rabani, E.; Baer, R. Distribution of Multiexciton Generation Rates in CdSe and InAs Nanocrystals. Nano Lett. 2008, 8, 4488-4492. 189. Kilina, S. V.; Kilin, D. S.; Prezhdo, O. V. Breaking the Phonon Bottleneck in PbSe and CdSe Quantum Dots: Time-Domain Density Functional Theory of Charge Carrier Relaxation. ACS Nano 2009, 3, 93-99. 190. Li, L. Q.; Long, R.; Prezhdo, O. V. Charge Separation and Recombination in Two- Dimensional MoS2/WS2: Time-Domain ab Initio Modeling. Chem. Mater. 2017, 29, 2466-2473. 191. Chu, W. B.; Saidi, W. A.; Prezhdo, O. V. Long-Lived Hot Electron in a Metallic Particle for Plasmonics and Catalysis: Ab Initio Nonadiabatic Molecular Dynamics with Machine Learning. ACS Nano 2020, 14, 10608-10615. 192. Li, W.; Vasenko, A. S.; Tang, J. F.; Prezhdo, O. V. Anharmonicity Extends Carrier Lifetimes in Lead Halide Perovskites at Elevated Temperatures. J. Phys. Chem. Lett. 2019, 10, 6219-6226. 193. Rong, Y. G.; Hu, Y.; Mei, A. Y.; Tan, H. R.; Saidaminov, M. I.; Seok, S. I.; McGehee, M. D.; Sargent, E. H.; Han, H. W. Challenges for commercializing perovskite solar cells. Science 2018, 361. 194. Xing, G.; Mathews, N.; Sun, S.; Lim, S. S.; Lam, Y. M.; Grätzel, M.; Mhaisalkar, S.; Sum, T. C. Long-range balanced electron-and hole-transport lengths in organic-inorganic CH3NH3PbI3. Science 2013, 342, 344-347. 195. Green, M. A.; Ho-Baillie, A. Perovskite solar cells: the birth of a new era in photovoltaics. ACS Energy Lett. 2017, 2, 822-830. 196. Manser, J. S.; Christians, J. A.; Kamat, P. V. Intriguing optoelectronic properties of metal halide perovskites. Chem. Rev. 2016, 116, 12956-13008. 197. Brenner, T. M.; Egger, D. A.; Kronik, L.; Hodes, G.; Cahen, D. Hybrid organic—inorganic perovskites: low-cost semiconductors with intriguing charge-transport properties. Nat. Rev. Mater. 2016, 1, 15007. 198. Herz, L. M. Charge-Carrier Mobilities in Metal Halide Perovskites: Fundamental Mechanisms and Limits. ACS Energy Lett. 2017, 2, 1539-1548. 133 199. Yin, W. J.; Yang, J. H.; Kang, J.; Yan, Y. F.; Wei, S. H. Halide perovskite materials for solar cells: a theoretical review. J. Mat. Chem. A 2015, 3, 8926-8942. 200. Liang, P. W.; Liao, C. Y.; Chueh, C. C.; Zuo, F.; Williams, S. T.; Xin, X. K.; Lin, J. J.; Jen, A. K. Y. Additive Enhanced Crystallization of Solution-Processed Perovskite for Highly Efficient Planar-Heterojunction Solar Cells. Adv. Mater. 2014, 26, 3748-3754. 201. Wu, Y.; Li, X.; Zeng, H. Highly luminescent and stable halide perovskite nanocrystals. ACS Energy Lett. 2019, 4, 673-681. 202. Yang, W. S.; Park, B.-W.; Jung, E. H.; Jeon, N. J.; Kim, Y. C.; Lee, D. U.; Shin, S. S.; Seo, J.; Kim, E. K.; Noh, J. H. Iodide management in formamidinium-lead-halide–based perovskite layers for efficient solar cells. Science 2017, 356, 1376-1379. 203. Tan, Z. K.; Moghaddam, R. S.; Lai, M. L.; Docampo, P.; Higler, R.; Deschler, F.; Price, M.; Sadhanala, A.; Pazos, L. M.; Credgington, D.; Hanusch, F.; Bein, T.; Snaith, H. J.; Friend, R. H. Bright light-emitting diodes based on organometal halide perovskite. Nat. Nanotechnol. 2014, 9, 687-692. 204. Li, G. R.; Tan, Z. K.; Di, D. W.; Lai, M. L.; Jiang, L.; Lim, J. H. W.; Friend, R. H.; Greenham, N. C. Efficient Light-Emitting Diodes Based on Nanocrystalline Perovskite in a Dielectric Polymer Matrix. Nano Lett. 2015, 15, 2640-2644. 205. Yantara, N.; Bhaumik, S.; Yan, F.; Sabba, D.; Dewi, H. A.; Mathews, N.; Boix, P. P.; Demir, H. V.; Mhaisalkar, S. Inorganic Halide Perovskites for Efficient Light-Emitting Diodes. J. Phys. Chem. Lett. 2015, 6, 4360-4364. 206. Miyano, K.; Tripathi, N.; Yanagida, M.; Shirai, Y. Lead Halide Perovskite Photovoltaic as a Model p-i-n Diode. Acc. Chem. Res. 2016, 49, 303-310. 207. Huang, J.; Yuan, Y.; Shao, Y.; Yan, Y. Understanding the physical properties of hybrid perovskites for photovoltaic applications. Nat. Rev. Mater. 2017, 2. 208. Luo, D. Y.; Yang, W. Q.; Wang, Z. P.; Sadhanala, A.; Hu, Q.; Su, R.; Shivanna, R.; Trindade, G. F.; Watts, J. F.; Xu, Z. J.; Liu, T. H.; Chen, K.; Ye, F. J.; Wu, P.; Zhao, L. C.; Wu, J.; Tu, Y. G.; Zhang, Y. F.; Yang, X. Y.; Zhang, W.; Friend, R. H.; Gong, Q. H.; Snaith, H. J.; Zhu, R. Enhanced photovoltage for inverted planar heterojunction perovskite solar cells. Science 2018, 360, 1442-1446. 209. Parvazian, E.; Abdollah-Zadeh, A.; Akbari, H. R.; Taghavinia, N. Fabrication of perovskite solar cells based on vacuum-assisted linear meniscus printing of MAPbI(3). Sol. Energy Mater. Sol. Cells 2019, 191, 148-156. 210. Correa-Baena, J.-P.; Saliba, M.; Buonassisi, T.; Grätzel, M.; Abate, A.; Tress, W.; Hagfeldt, A. Promises and challenges of perovskite solar cells. Science 2017, 358, 739-744. 211. Kwon, H. C.; Yang, W.; Lee, D.; Ahn, J.; Lee, E.; Ma, S.; Kim, K.; Yun, S. C.; Moon, J. Investigating Recombination and Charge Carrier Dynamics in a One-Dimensional Nanopillared Perovskite Absorber. ACS Nano 2018, 12, 4233-4245. 212. Fang, Y.; Bi, C.; Wang, D.; Huang, J. The functions of fullerenes in hybrid perovskite solar cells. ACS Energy Lett. 2017, 2, 782-794. 213. Li, Z.; Kolodziej, C.; McCleese, C.; Wang, L.; Kovalsky, A.; Samia, A. C.; Zhao, Y.; Burda, C. Effect of chloride substitution on interfacial charge transfer processes in MAPbI3 perovskite thin film solar cells: planar versus mesoporous. Nanoscale Adv. 2019, 1, 827-833. 214. Aranda, C.; Guerrero, A.; Bisquert, J. Ionic Effect Enhances Light Emission and Photovoltage of Methylammonium Lead Bromide Perovskite Solar Cell by Reduced Surface Recombination. ACS Energy Lett. 2019. 134 215. Liu, M. Z.; Johnston, M. B.; Snaith, H. J. Efficient planar heterojunction perovskite solar cells by vapour deposition. Nature 2013, 501, 395-398. 216. Zhou, H. P.; Chen, Q.; Li, G.; Luo, S.; Song, T. B.; Duan, H. S.; Hong, Z. R.; You, J. B.; Liu, Y. S.; Yang, Y. Interface engineering of highly efficient perovskite solar cells. Science 2014, 345, 542-546. 217. Yang, W. S.; Noh, J. H.; Jeon, N. J.; Kim, Y. C.; Ryu, S.; Seo, J.; Seok, S. I. High- performance photovoltaic perovskite layers fabricated through intramolecular exchange. Science 2015, 348, 1234-1237. 218. Burschka, J.; Pellet, N.; Moon, S. J.; Humphry-Baker, R.; Gao, P.; Nazeeruddin, M. K.; Gratzel, M. Sequential deposition as a route to high-performance perovskite-sensitized solar cells. Nature 2013, 499, 316-319. 219. Du, M. H. Efficient carrier transport in halide perovskites: theoretical perspectives. J. Mat. Chem. A 2014, 2, 9091-9098. 220. He, J. L.; Fang, W. H.; Long, R.; Prezhdo, O. V. Superoxide/Peroxide Chemistry Extends Charge Carriers' Lifetime but Undermines Chemical Stability of CH3NH3PbI3 Exposed to Oxygen: Time-Domain ab Initio Analysis. J. Am. Chem. Soc. 2019, 141, 5798-5807. 221. Chakraborty, S.; Xie, W.; Mathews, N.; Sherburne, M.; Ahuja, R.; Asta, M.; Mhaisalkar, S. G. Rational Design: A High-Throughput Computational Screening and Experimental Validation Methodology for Lead-Free and Emergent Hybrid Perovskites. ACS Energy Lett. 2017, 2, 837- 845. 222. Lo, Y. C.; Rensi, S. E.; Torng, W.; Altman, R. B. Machine learning in chemoinformatics and drug discovery. Drug Discovery Today 2018, 23, 1538-1546. 223. Begone, G.; Deisenroth, M. P.; Kim, J. S.; Liem, S.; de Austri, R. R.; Welling, M. Accelerating the BSM interpretation of LHC data with machine learning. Phys. Dark Universe 2019, 24. 224. Yu, Y. Z.; Tan, X. H.; Ning, S. G.; Wu, Y. Y. Machine Learning for Understanding Compatibility of Organic-Inorganic Hybrid Perovskites with Post-Treatment Amines. ACS Energy Lett. 2019, 4, 397-404. 225. Butler, K. T.; Davies, D. W.; Cartwright, H.; Isayev, O.; Walsh, A. Machine learning for molecular and materials science. Nature 2018, 559, 547-555. 226. Ghasemi, S. A.; Hofstetter, A.; Saha, S.; Goedecker, S. Interatomic potentials for ionic systems with density functional accuracy based on charge densities obtained by a neural network. Phys. Rev. B 2015, 92, 045131. 227. Tavadze, P.; Avendañ o Franco, G.; Ren, P.; Wen, X.; Li, Y.; Lewis, J. P. A machine-driven hunt for global reaction coordinates of azobenzene photoisomerization. J. Am. Chem. Soc. 2017, 140, 285-290. 228. Prezhdo, O. V. Multiple excitons and the electron-phonon bottleneck in semiconductor quantum dots: An ab initio perspective. Chem. Phys. Lett. 2008, 460, 1-9. 229. Zhou, X.; Jankowska, J.; Dong, H.; Prezhdo, O. V. Recent theoretical progress in the development of perovskite photovoltaic materials. J. Ener. Chem. 2018, 27, 637-649. 230. Long, R.; Prezhdo, O. V.; Fang, W. H. Nonadiabatic charge dynamics in novel solar cell materials. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2017, 7. 231. Li, W.; Tang, J. F.; Casanova, D.; Prezhdo, O. V. Time-Domain ab lnitio Analysis Rationalizes the Unusual Temperature Dependence of Charge Carrier Relaxation in Lead Halide Perovskite. ACS Energy Lett. 2018, 3, 2713-2720. 135 232. Leguy, A. M.; Azarhoosh, P.; Alonso, M. I.; Campoy-Quiles, M.; Weber, O. J.; Yao, J.; Bryant, D.; Weller, M. T.; Nelson, J.; Walsh, A. Experimental and theoretical optical properties of methylammonium lead halide perovskites. Nanoscale 2016, 8, 6317-6327. 233. Baikie, T.; Fang, Y.; Kadro, J. M.; Schreyer, M.; Wei, F.; Mhaisalkar, S. G.; Graetzel, M.; White, T. J. Synthesis and crystal chemistry of the hybrid perovskite (CH 3 NH 3) PbI 3 for solid- state sensitised solar cell applications. J. Mat. Chem. A 2013, 1, 5628-5641. 234. Zhou, X.; Jankowska, J.; Li, L. Q.; Giri, A.; Hopkins, P. E.; Prezhdo, O. V. Strong Influence of Ti Adhesion Layer on Electron-Phonon Relaxation in Thin Gold Films: Ab Initio Nonadiabatic Molecular Dynamics. ACS Appl. Mater. Interfaces 2017, 9, 43343-43351. 235. Akimov, A. V.; Muckerman, J. T.; Prezhdo, O. V. Nonadiabatic Dynamics of Positive Charge during Photocatalytic Water Splitting on GaN(10-10) Surface: Charge Localization Governs Splitting Efficiency. J. Am. Chem. Soc. 2013, 135, 8682-8691. 236. Chaban, V. V.; Prezhdo, V. V.; Prezhdo, O. V. Covalent Linking Greatly Enhances Photoinduced Electron Transfer in Fullerene-Quantum Dot Nanocomposites: Time-Domain Ab Initio Study. J. Phys. Chem. Lett. 2013, 4, 1-6. 237. He, J.; Fang, W.-H.; Long, R.; Prezhdo, O. V. Increased Lattice Stiffness Suppresses Nonradiative Charge Recombination in MAPbI3 Doped with Larger Cations: Time-Domain Ab Initio Analysis. ACS Energy Lett. 2018, 3, 2070-2076. 238. Li, W.; Liu, J.; Bai, F.-Q.; Zhang, H.-X.; Prezhdo, O. V. Hole trapping by iodine interstitial defects decreases free carrier losses in perovskite solar cells: a time-domain ab initio study. ACS Energy Lett. 2017, 2, 1270-1278. 239. Kraskov, A.; Stogbauer, H.; Grassberger, P. Estimating mutual information. Phys. Rev. E 2004, 69. 240. Khan, S.; Bandyopadhyay, S.; Ganguly, A. R.; Saigal, S.; Erickson, D. J.; Protopopescu, V.; Ostrouchov, G. Relative performance of mutual information estimation methods for quantifying the dependence among short and noisy data. Phys. Rev. E 2007, 76. 241. Hyeon-Deuk, K.; Madrid, A. B.; Prezhdo, O. V. Symmetric band structures and asymmetric ultrafast electron and hole relaxations in silicon and germanium quantum dots: time- domain ab initio simulation. Dalton Trans. 2009, 10069-10077. 242. Mahata, A.; Meggiolaro, D.; De Angelis, F. From Large to Small Polarons in Lead, Tin, and Mixed Lead-Tin Halide Perovskites. J. Phys. Chem. Lett. 2019, 10, 1790-1798. 243. Miyata, K.; Meggiolaro, D.; Trinh, M. T.; Joshi, P. P.; Mosconi, E.; Jones, S. C.; De Angelis, F.; Zhu, X. Y. Large polarons in lead halide perovskites. Science Adv. 2017, 3. 244. Milot, R. L.; Eperon, G. E.; Snaith, H. J.; Johnston, M. B.; Herz, L. M. Temperature- Dependent Charge-Carrier Dynamics in CH3NH3PbI3 Perovskite Thin Films. Adv. Funct. Mater. 2015, 25, 6218-6227. 245. Ambrosio, F.; Wiktor, J.; De Angelis, F.; Pasquarello, A. Origin of low electron-hole recombination rate in metal halide perovskites. Energy Environ. Sci. 2018, 11, 101-105. 246. Neukirch, A. J.; Nie, W. Y.; Blancon, J. C.; Appavoo, K.; Tsai, H.; Sfeir, M. Y.; Katan, C.; Pedesseau, L.; Even, J.; Crochet, J. J.; Gupta, G.; Mohite, A. D.; Tretiak, S. Polaron Stabilization by Cooperative Lattice Distortion and Cation Rotations in Hybrid Perovskite Materials. Nano Lett. 2016, 16, 3809-3816. 247. Park, M.; Neukirch, A. J.; Reyes-Lillo, S. E.; Lai, M. L.; Ellis, S. R.; Dietze, D.; Neaton, J. B.; Yang, P. D.; Tretiak, S.; Mathies, R. A. Excited-state vibrational dynamics toward the polaron in methylammonium lead iodide perovskite. Nat. Commun. 2018, 9. 136 248. Even, J.; Pedesseau, L.; Katan, C. Analysis of Multivalley and Multibandgap Absorption and Enhancement of Free Carriers Related to Exciton Screening in Hybrid Perovskites. J. Phys. Chem. C 2014, 118, 11566-11572. 249. Zhu, H. M.; Miyata, K.; Fu, Y. P.; Wang, J.; Joshi, P. P.; Niesner, D.; Williams, K. W.; Jin, S.; Zhu, X. Y. Screening in crystalline liquids protects energetic carriers in hybrid perovskites. Science 2016, 353, 1409-1413. 250. Tong, C. J.; Geng, W.; Prezhdo, O. V.; Liu, L. M. Role of Methylammonium Orientation in Ion Diffusion and Current Voltage Hysteresis in the CH3NH3PbI3 Perovskite. ACS Energy Lett. 2017, 2, 1997-2004. 251. Lahnsteiner, J.; Kresse, G.; Kumar, A.; Sarma, D.; Franchini, C.; Bokdam, M. Room- temperature dynamic correlation between methylammonium molecules in lead-iodine based perovskites: An ab initio molecular dynamics perspective. Phys. Rev. B 2016, 94, 214114. 252. Thiel, W. Semiempirical quantum–chemical methods. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2014, 4, 145-157. 253. Bikadi, Z.; Hazai, E. Application of the PM6 semi-empirical method to modeling proteins enhances docking accuracy of AutoDock. J. Cheminf. 2009, 1, 15. 254. Dewar, M. J.; Thiel, W. Ground states of molecules. 38. The MNDO method. Approximations and parameters. J. Am. Chem. Soc. 1977, 99, 4899-4907. 255. Dewar, M. J.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. Development and use of quantum mechanical molecular models. 76. AM1: a new general purpose quantum mechanical molecular model. J. Am. Chem. Soc. 1985, 107, 3902-3909. 256. Stewart, J. J. Optimization of parameters for semiempirical methods II. Applications. J. Comput. Chem. 1989, 10, 221-264. 257. Stewart, J. J. Optimization of parameters for semiempirical methods V: modification of NDDO approximations and application to 70 elements. J. Mol. Model. 2007, 13, 1173-1213. 258. Dral, P. O.; Wu, X.; Thiel, W. Semiempirical quantum-chemical methods with Orthogonalization and dispersion corrections. J. Chem. Theory Comput. 2019, 15, 1743-1760. 259. McNamara, J. P.; Sharma, R.; Vincent, M. A.; Hillier, I. H.; Morgado, C. A. The non- covalent functionalisation of carbon nanotubes studied by density functional and semi-empirical molecular orbital methods including dispersion corrections. Phys. Chem. Chem. Phys. 2008, 10, 128-135. 260. Sure, R.; Grimme, S. Comprehensive benchmark of association (free) energies of realistic host–guest complexes. J. Chem. Theory Comput. 2015, 11, 3785-3801. 261. Dobeš , P.; Rezác, J.; Fanfrlik, J.; Otyepka, M.; Hobza, P. Semiempirical quantum mechanical method PM6-DH2X describes the geometry and energetics of CK2-inhibitor complexes involving halogen bonds well, while the empirical potential fails. J. Phys. Chem. B 2011, 115, 8581-8589. 262. Thiel, W.; Voityuk, A. A. Extension of MNDO to d orbitals: parameters and results for silicon. J. Mol. Struct.: THEOCHEM 1994, 313, 141-154. 263. Stewart, J. J. Optimization of parameters for semiempirical methods VI: more modifications to the NDDO approximations and re-optimization of parameters. J. Mol. Model. 2013, 19, 1-32. 264. Rezac, J.; Fanfrlik, J.; Salahub, D.; Hobza, P. Semiempirical quantum chemical PM6 method augmented by dispersion and H-bonding correction terms reliably describes various types of noncovalent complexes. J. Chem. Theory Comput. 2009, 5, 1749-1760. 137 265. Korth, M.; Pitonak, M.; Rezac, J.; Hobza, P. A transferable H-bonding correction for semiempirical quantum-chemical methods. J. Chem. Theory Comput. 2009, 6, 344-352. 266. Korth, M. Third-generation hydrogen-bonding corrections for semiempirical QM methods and force fields. J. Chem. Theory Comput. 2010, 6, 3808-3816. 267. Dral, P. O.; Wu, X.; Spö rkel, L.; Koslowski, A.; Weber, W.; Steiger, R.; Scholten, M.; Thiel, W. Semiempirical quantum-chemical orthogonalization-corrected methods: theory, implementation, and parameters. J. Chem. Theory Comput. 2016, 12, 1082-1096. 268. Grimme, S.; Steinmetz, M. Effects of London dispersion correction in density functional theory on the structures of organic molecules in the gas phase. Phys. Chem. Chem. Phys. 2013, 15, 16031-16042. 269. Bereau, T.; DiStasio, R. A.; Tkatchenko, A.; von Lilienfeld, O. A. Non-covalent interactions across organic and biological subsets of chemical space: Physics-based potentials parametrized from machine learning. J. Chem. Phys. 2018, 148, 241706. 270. Niklasson, A. M. Expansion algorithm for the density matrix. Phys. Rev. B 2002, 66, 155115. 271. Paszke, A.; Gross, S.; Chintala, S.; Chanan, G.; Yang, E.; DeVito, Z.; Lin, Z.; Desmaison, A.; Antiga, L.; Lerer, A. Automatic differentiation in pytorch. NIPS 2017 Workshop Autodif 2017. 272. Ahneman, D. T.; Estrada, J. G.; Lin, S.; Dreher, S. D.; Doyle, A. G. Predicting reaction performance in C–N cross-coupling using machine learning. Science 2018, 360, 186-190. 273. Raccuglia, P.; Elbert, K. C.; Adler, P. D.; Falk, C.; Wenny, M. B.; Mollo, A.; Zeller, M.; Friedler, S. A.; Schrier, J.; Norquist, A. J. Machine-learning-assisted materials discovery using failed experiments. Nature 2016, 533, 73. 274. Lubbers, N.; Smith, J. S.; Barros, K. Hierarchical modeling of molecular energies using a deep neural network. J. Chem. Phys. 2018, 148, 241715. 275. Zubatyuk, T.; Nebgen, B.; Lubbers, N.; Smith, J. S.; Zubatyuk, R.; Zhou, G.; Koh, C.; Barros, K.; Isayev, O.; Tretiak, S. Machine Learned Hückel Theory: Interfacing Physics and Deep Neural Networks. arXiv preprint:1909.12963 2019. 276. Zhou, G. PYSEQM. https://github.com/lanl/PYSEQM (accessed June 2020). 277. Tapavicza, E.; Bellchambers, G. D.; Vincent, J. C.; Furche, F. Ab initio non-adiabatic molecular dynamics. Phys. Chem. Chem. Phys. 2013, 15, 18336-18348. 278. Niklasson, A. M. Extended born-oppenheimer molecular dynamics. Phys. Rev. Lett. 2008, 100, 123004. 279. Niklasson, A. M.; Steneteg, P.; Odell, A.; Bock, N.; Challacombe, M.; Tymczak, C.; Holmström, E.; Zheng, G.; Weber, V. Extended Lagrangian Born–Oppenheimer molecular dynamics with dissipation. J. Chem. Phys. 2009, 130, 214109. 280. Souvatzis, P.; Niklasson, A. M. First principles molecular dynamics without self-consistent field optimization. J. Chem. Phys. 2014, 140, 044117. 281. Niklasson, A. M.; Cawkwell, M. J. Generalized extended lagrangian born-oppenheimer molecular dynamics. J. Chem. Phys. 2014, 141, 164123. 282. Niklasson, A. M. Next generation extended Lagrangian first principles molecular dynamics. J. Chem. Phys. 2017, 147, 054103. 283. Mniszewski, S. M.; Cawkwell, M. J.; Wall, M. E.; Mohd-Yusof, J.; Bock, N.; Germann, T. C.; Niklasson, A. Efficient parallel linear scaling construction of the density matrix for Born– Oppenheimer molecular dynamics. J. Chem. Theory Comput. 2015, 11, 4644-4654. 138 284. Cawkwell, M.; Wood, M.; Niklasson, A. M.; Mniszewski, S. Computation of the density matrix in electronic structure theory in parallel on multiple graphics processing units. J. Chem. Theory Comput. 2014, 10, 5391-5396. 285. Rubensson, E. H.; Niklasson, A. M. Interior eigenvalues from density matrix expansions in quantum mechanical molecular dynamics. SIAM J. Sci. Comput. 2014, 36, B147-B170. 286. Mniszewski, S. M.; Perriot, R.; Rubensson, E. H.; Negre, C. F. A.; Cawkwell, M. J.; Niklasson, A. M. N. Linear Scaling Pseudo Fermi-Operator Expansion for Fractional Occupation. J. Chem. Theory Comput. 2018, 15, 190-200. 287. Palma, J. L.; Atas, E.; Hardison, L.; Marder, T. B.; Collings, J. C.; Beeby, A.; Melinger, J. S.; Krause, J. L.; Kleiman, V. D.; Roitberg, A. E. Electronic spectra of the nanostar dendrimer: Theory and experiment. J. Phys. Chem. C 2010, 114, 20702-20712. 288. Roothaan, C. C. J. New developments in molecular orbital theory. Rev. Mod. Phys. 1951, 23, 69. 289. Dewar, M. J.; Thiel, W. A semiempirical model for the two-center repulsion integrals in the NDDO approximation. Theor. Chim. Acta 1977, 46, 89-104. 290. Badziag, P.; Solms, F. An improved SCF iteration scheme. Comput. Chem. 1988, 12, 233- 236. 291. Pulay, P. Convergence acceleration of iterative sequences. The case of SCF iteration. Chem. Phys. Lett. 1980, 73, 393-398. 292. Camp, R. N.; King, H. F. An interpolation method for forcing SCF convergence. J. Chem. Phys. 1981, 75, 268-274. 293. Chaban, G.; Schmidt, M. W.; Gordon, M. S. Approximate second order method for orbital optimization of SCF and MCSCF wavefunctions. Theor. Chem. Acc. 1997, 97, 88-95. 294. Daniels, A. D.; Scuseria, G. E. What is the best alternative to diagonalization of the Hamiltonian in large scale semiempirical calculations? J. Chem. Phys. 1999, 110, 1321-1328. 295. Niklasson, A. M.; Tymczak, C.; Challacombe, M. Time-reversible Born-Oppenheimer molecular dynamics. Phys. Rev. Lett. 2006, 97, 123001. 296. Verlet, L. Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Phys. Rev. 1967, 159, 98. 297. Neese, F. The ORCA program system. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2012, 2, 73-78. 298. Fernandez-Alberti, S.; Kleiman, V. D.; Tretiak, S.; Roitberg, A. E. Nonadiabatic molecular dynamics simulations of the energy transfer between building blocks in a phenylene ethynylene dendrimer. J. Phys. Chem. A 2009, 113, 7535-7542. 299. Bradshaw, D. S.; Andrews, D. L. Mechanisms of light energy harvesting in dendrimers and hyperbranched polymers. Polymers 2011, 3, 2053-2077. 300. Stuart, S. J.; Tutein, A. B.; Harrison, J. A. A reactive potential for hydrocarbons with intermolecular interactions. J. Chem. Phys. 2000, 112, 6472-6486. 301. Neese, F. Approximate second-order SCF convergence for spin unrestricted wavefunctions. Chem. Phys. Lett. 2000, 325, 93-98. 302. Hamilton, T. P.; Pulay, P. Direct inversion in the iterative subspace (DIIS) optimization of open ‐shell, excited ‐state, and small multiconfiguration SCF wave functions. J. Chem. Phys. 1986, 84, 5728-5734. 303. Toxvaerd, S.; Heilmann, O. J.; Dyre, J. C. Energy conservation in molecular dynamics simulations of classical systems. J. Chem. Phys. 2012, 136, 224106. 139 304. Lippert, R. A.; Bowers, K. J.; Dror, R. O.; Eastwood, M. P.; Gregersen, B. A.; Klepeis, J. L.; Kolossvary, I.; Shaw, D. E. A common, avoidable source of error in molecular dynamics integrators. J. Chem. Phys. 2007, 126, 046101. 305. Zhou, G.; Cen, C.; Wang, S.; Deng, M.; Prezhdo, O. V. Electron-Phonon Scattering is Much Weaker in Carbon Nanotubes than Graphene Nanoribbons. J. Phys. Chem. Lett. 2019. 306. Callis, P. R.; Vivian, J. T. Understanding the variable fluorescence quantum yield of tryptophan in proteins using QM-MM simulations. Quenching by charge transfer to the peptide backbone. Chem. Phys. Lett. 2003, 369, 409-414. 307. Shalashilin, D. V. Nonadiabatic dynamics with the help of multiconfigurational Ehrenfest method: Improved theory and fully quantum 24D simulation of pyrazine. J. Chem. Phys. 2010, 132, 244111. 308. Ramezanghorbani, F.; Isayev, O.; Smith, J.; Roitberg, A. TorchANI: A Free and Open Source PyTorch Based Deep Learning Implementation of the ANI Neural Network Potentials. 2020. 309. Schütt, K. T.; Sauceda, H. E.; Kindermans, P.-J.; Tkatchenko, A.; Müller, K.-R. SchNet– A deep learning architecture for molecules and materials. J. Chem. Phys. 2018, 148, 241722. 310. Schütt, K.; Kindermans, P.-J.; Felix, H. E. S.; Chmiela, S.; Tkatchenko, A.; Müller, K.-R. Schnet: A continuous-filter convolutional neural network for modeling quantum interactions, Adv. Neural Inf. Process. Syst. , 2017; 991-1001. 311. Schutt, K.; Kessel, P.; Gastegger, M.; Nicoli, K.; Tkatchenko, A.; Muller, K.-R. SchNetPack: A deep learning toolbox for atomistic systems. J. Chem. Theory Comput. 2018, 15, 448-455. 312. Qiao, Z.; Welborn, M.; Anandkumar, A.; Manby, F. R.; Miller III, T. F. OrbNet: Deep learning for quantum chemistry using symmetry-adapted atomic-orbital features. J. Chem. Phys. 2020, 153, 124111. 313. Li, H.; Collins, C.; Tanha, M.; Gordon, G. J.; Yaron, D. J. A density functional tight binding layer for deep learning of chemical Hamiltonians. J. Chem. Theory Comput. 2018, 14, 5764-5776. 314. Smith, J. S.; Nebgen, B.; Lubbers, N.; Isayev, O.; Roitberg, A. E. Less is more: Sampling chemical space with active learning. J. Chem. Phys. 2018, 148, 241733. 315. Zhou, G.; Nebgen, B.; Lubbers, N.; Malone, W.; Niklasson, A. M.; Tretiak, S. Graphics processing unit-accelerated semiempirical Born Oppenheimer molecular dynamics using PyTorch. J. Chem. Theory Comput. 2020, 16, 4951-4962. 316. Frisch, M.; Trucks, G.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. gaussian 09, Revision d. 01, Gaussian. Inc., Wallingford CT 2009, 201. 317. Malone, W.; Nebgen, B.; White, A.; Zhang, Y.; Song, H.; Bjorgaard, J. A.; Sifain, A. E.; Rodriguez-Hernandez, B.; Freixas, V. M.; Fernandez-Alberti, S. NEXMD Software Package for Nonadiabatic Excited State Molecular Dynamics Simulations. J. Chem. Theory Comput. 2020, 16, 5771-5783. 318. Yosipof, A.; Nahum, O. E.; Anderson, A. Y.; Barad, H. N.; Zaban, A.; Senderowitz, H. Data Mining and Machine Learning Tools for Combinatorial Material Science of All ‐Oxide Photovoltaic Cells. Mol. Inf. 2015, 34, 367-379. 140 Appendix A: Supplementary Information for Chapter 2 A1. Force-Field Validation by Experimental Measurement of Contact Angle The simulation reported here was performed with a combination of force fields: REBO (reactive bond order) for MoS2 and TIP4P/2005 for H2O. REBO potential accounts for changes in local atomic configurations of Mo and S, and TIP4P/2005 correctly models structural and dynamical properties of water. The interaction between MoS2 and water is described by a combination of Lennard-Jones (LJ) and electrostatic potentials with force-field parameters taken from Luan et al. 61 We apply the Lorentz-Berthelot combination rule to parameterize LJ interactions between H2O, MoS2 and IPA molecules. Electrostatic forces and energy are calculated with the PPPM method, and slab correction 62 is applied to allow for fixed boundaries in the z direction. We use the Velocity-Verlet integrator with constraints on O-H bond and H-O-H bond angle of water molecules imposed with the SHAKE algorithm 19 . Force fields were validated by experimental data on contact angles and surface tensions 46, 63 . According to Young ’s law, this ensures the correctness of interfacial energy between the solvent and MoS2. The validation simulation was performed for a system consisting of a droplet of H2O/IPA mixture resting on an MoS2 substrate, see Figure S2.1. We use free boundary conditions along x and z and periodic boundary condition in the y direction. The droplet of solvent molecules was put on top of a four-layer thick MoS2 substrate and the whole system was relaxed under ambient conditions until the contact angle had converged. The initial force field parameters between IPA and MoS2 were generated using the mixing rule. The energy parameter  in the LJ force field was scaled by a parameter to fit two experimental measurements of the contact angle: one without IPA and the other for mass concentration of IPA around 50%. The final results are shown in the Table S2.1. 141 Mass fraction MD () Experiment () 0.0 1014 97.8 0.035 753 70.6 0.16 392 44.8 0.36 314 31.1 0.56 191 24.8 Table S2.1: Contact angles calculated by MD simulations for various mass fractions of IPA in the mixture. Experimental data are taken from Halim, U. et al 63 . After relaxing the system for 3 ns, six different configurations separated by 2 ps were used to calculate contact angles. Figure S2.1: Side view of the system for the calculation of contact angle of H2O/IPA (C, O: red, H: white) mixture on top of an MoS2 substrate (Mo: pink, S: yellow). Here the mass concentration of IPA is 36%. The calculated contact angle is 42 on the left and 36 on the right-hand side. A2. Force-Field Validation for Shock Simulation To get the shock hugoniot of the H2O/IPA mixture, we apply different particle velocities Vp in our H2O/IPA mixture system to get the shock velocity Vs. The system’s dimensions are 19.7nm×19.7nm×29.6nm in the x, y, and z directions, it contains 1.17 million atoms with 177332 water and 52974 IPA molecules (1:1 ratio by weight). The system is periodic along x and y directions while fixed along z with a momentum mirror on the ends of z direction. The system is first relaxed for 1 ns with timestep 2 fs under ambient condition. After that, a constant particle 142 velocity Vp is added to the whole system. When atoms crossing the momentum mirror, they are reflected with their momenta in the z direction reversed, which generate a planer shock wave in the mixture propagating into the system. The shock velocity Vs is calculated with shock front boundaries at different time frames. The boundary is located with the abrupt of density. The simulation results show that the water/IPA mixture has similar shock hugoniot as pure water, see Figure S2.2. Figure S2.2 shows the Hugoniot compression curves of water (red circles: MD simulation, blue crosses: experiment) and water/IPA mixture (green segments: MD simulation) (50 wt% of IPA). The MD results and the experimental data for pure water are from Vedadi, M. et al 52 A3. Additional Results To quantify the exfoliation yield during the simulation, the accessible surface area 64, 66 and the volume of convex hull 67 of the MoS2 bulk are calculated during the simulation as shown in Figure 143 S2.3(b). The convex hull of the bulk is used to count the number of solvent molecules flowing into and out of the MoS2 sheets during the shock wave propagation (see Figure S2.3(a) and (c)). Each time when the shock wave hits (pass through) the bulk, solvent molecules are flowing into (out of) the sheets. Figure S2.3: Exfoliated MoS2 with solvent between the nanosheets. (a) Exfoliated MoS2 with only a fraction (10%) of solvent molecules shown in the figure at t = 40ps. (b) shows changes in the volume (blue line) and surface area (red line) of MoS2 during exfoliation as a function of time. The volume is that of the convex hull computed from the largest clusters of MoS2. The surface area is computed from the surface mesh of the bulk. Both the volume and surface area are normalized by their initial values. At the end of exfoliation the MoS2 volume expands by 100%, and the surface area increases by 60%. (c) shows the percentage of water (red line) and IPA (blue line) molecules in the convex hull as a function of time. Each time the shock wave hits MoS2 (at t = 2 ps, 12 ps, 18 ps), the solvent flows into the galleries of MoS2. After the shock wave passes, some of the solvent molecules flow out of the galleries. Our simulations show that the exfoliation happens for standoff parameter S between 1.1 to 2.0. In Figure S2.4 (a) we show exfoliation of MoS2 for S = 2.0. We performed simulations for particle velocities Vp = 0.5, 1.0, 2.0, 3.0, 4.0 km/s and did not observe exfoliation for Vp < 3 km/s. At Vp = 4.0 km/s, the nanojet impact disorders MoS2 and increases its temperature to 3,500 K. This indicates that the intensity of sonication plays a crucial role in exfoliation. We also performed simulations at two different orientations of MoS2, one in which <0001> and the other in which <21 ̅ 1 ̅ 0> surface of MoS2 is oriented to normal to shock propagation. At Vp 144 = 3.0 km/s we find that the nanojet impact on the <0001> surface yields better exfoliation than the impact on the <21 ̅ 1 ̅ 0> surface (see Figure S2.4 (b)). Figure S2.4 shows exfoliation of MoS2 samples (Mo: pink, S: yellow) in two MD simulations. (a) shows exfoliated MoS2 for S = 2.0, bubble size 5.0 nm, and particle velocity Vp 3.0 km/s when the shock is normal to the <21 ̅ 1 ̅ 0> surface. (b) gives side view of cleavage in bulk MoS2 when the nanojet impacts the <0001> surface. In this case, the standoff parameter is 1.21, bubble size is 4.7 nm, and Vp is 3.0 km/s. 145 Appendix B: Supporting Information for Chapter 4 Hot electron and hole relaxation In the CdSe QD, the conduction band above the 1Pe state and the valence band near the band edge have dense energy spectra, and electron-phonon scattering dominates the carrier relaxation. Therefore, fewest switch surface hopping (FSSH) implemented in PYXAID is used to model hot electron and hole relaxation in these energy ranges. 10, 15 Figure S4.1 shows the decay of the energy of hot electrons initially excited above 1Pe. The curves are fitted with the exponential function Aexp(-t/), where  is the timescale parameter from fitting. The hot electron relaxation to the edge state 1Pe (LUMO+1) shows an average timescale of around 0.3 ps. Figure S4.2 presents the decay of the energy of hot holes excited below HOMO. The hole relaxation timescale is within the 0.2~0.3 ps range. Figure S4.1. Energy decay of electrons from states LUMO+(2,3,4,5) to 1Pe (LUMO+1). The reference zero energy is the energy of 1Pe. The timescales are obtained by exponential fits. 146 Figure S4.2. Energy decay of holes from states HOMO-(1, 2,…,14) to HOMO. The reference zero energy is the energy of HOMO. The timescales are obtained by exponential fits 147 Appendix C: Supporting Information for Chapter 5 C1.Theoretical Methodology In the simulation for the non-adiabatic molecular dynamics, we apply time-dependent density functional theory to describe the evolution of electrons 15 . The multielectron wavefunction Ψ(𝐫 ,𝐑 ,t) is governed by time-dependent Schrödinger equation (TDSE), 𝑖ℏ 𝜕 Ψ(𝐫 ,𝐑 ,t) 𝜕𝑡 =𝐻 (𝐫 ,𝑹 ,𝑡 )Ψ(𝐫 ,𝐑 ,t) (S5.1) Here H is the Hamiltonian describing the interaction between electrons and nuclei. Ψ is represented as a liner combination of basis functions Φ 𝑖 , which are slater determinants of adiabatic Kohn-Sham (KS) orbitals, Ψ(𝐫 ,𝐑 ,t)=∑𝑐 𝑖 (𝑡 )Φ 𝑖 (𝐫 |𝐑 ) 𝑖 (S5.2) The adiabatic KS orbitals can be obtained from any standard DFT packages. With this representation, the TDSE equation (S5.1) is expressed as the equations for expansion coefficients, 𝑖ℏ 𝜕 𝜕𝑡 𝑐 𝑖 (𝑡 )=∑(𝘀 𝑖 𝛿 𝑖𝑗 +𝑑 𝑖𝑗 )𝑐 𝑗 (𝑡 ) 𝑗 (S5.3) Where 𝘀 𝑖 is the energy associated with the basis function Φ 𝑖 , and 𝑑 𝑖𝑗 is the non-adiabatic coupling (NAC) between Φ 𝑖 and Φ 𝑗 , 𝑑 𝑖𝑗 = −𝑖ℏ<Φ 𝑖 | 𝜕 𝜕𝑡 |Φ 𝑗 > =− 𝑖ℏ<Φ 𝑖 |∇ 𝐑 |Φ 𝑗 > 𝐑 ̇ (S5.4) which can be reduced to non-adiabatic coupling between the KS orbitals. 𝑑 𝑖𝑗 is non-zero only if the configurations of Φ 𝑖 ,Φ 𝑗 differ by no more than one KS orbital. 15 𝑑 𝑖𝑗 is numerically computed as: 𝑑 𝑖𝑗 = −𝑖ℏ<Φ 𝑖 | 𝜕 𝜕𝑡 |Φ 𝑗 > ≈ − 𝑖 ℏ 2∆𝑡 (<Φ 𝑖 (𝑡 )|Φ 𝑗 (𝑡 +∆𝑡 )>−<Φ 𝑖 (𝑡 +∆𝑡 )|Φ 𝑗 (𝑡 )>) (S5.5) 148 C2. Features generated from the Molecular Dynamics trajectories. There are 22 static types of features and 22 corresponding dynamic ones included in the analysis for the mutual information computation with NAC, and the static ones are also used to the analysis with bandgap. These static features and the corresponding individual internal coordinates and motions are listed in Table S5.I. The dynamic features are the time derivative of the corresponding static features, which are computed for each individual coordinate as 𝑥 ̇(𝑡 + ∆𝑡 2 )=(𝑥 (𝑡 +∆𝑡 )−𝑥 (𝑡 ))/∆𝑡 (S5.6) Table S5.1: Static features used in the mutual information analysis Static Features Examples No. of individual ones in 111 system No. of individual ones in 222 system C-N bond C1-N1 in Fig. S5.1(a) 4 32 C-H bond C1-H1 in Fig. S5.1(a) 12 96 N-H bond N1-H13 in Fig. S5.1(a) 12 96 Pb-I bond Pb1-I2 in Fig. S5.1(b) 24 768 Pb-Pb distance 1 st and 2 nd nearest neighbor Pb1 -Pb2 in Fig. S5.1(b) 4 96 H-C-N bond angle H1-C1-N1 in Fig. S5.1(a) 4 96 H-C-H bond angle H1-C1-H2 in Fig. S5.1(a) 8 96 C-N-H bond angle C1-N1-H13 In Fig. S5.1(a) 4 96 H-N-H bond angle H13-N1-H14 in Fig. S5.1(a) 8 96 Pb-I-Pb bond angle Pb1-I3-Pb2 in Fig. S5.1(a) 12 96 I-I-I angle I6-I3-I5, I4-I3-I7, I10-I3-I9, I11-I3-I8 in Fig. S5.1(b) 15 256 I-Pb-I bond angle (flat) I2-Pb1-I3, I8-Pb1-I9, I5-Pb1-I7 in Fig. S5.1(b) 12 96 149 I-Pb-I bond angle (vertical) I2-Pb1-I9, I2-Pb1-I7, I2-Pb1-I5, I2-Pb1-I8 in Fig. S5.1(b) 48 384 H-C-N-H dihedral angle H1-C1-N1-H13 in Fig. S5. 1(a) 4 288 MA-axis c angle Angle  in Fig. S5.1(c) 4 32 MA-axis a angle Angle  in Fig. S5.1(c) 4 32 MA-MA angle Angle between C1- N1 and C3-N3 in Fig. S5.1(c) 6 256 MA MA center of mass distance Distance between P and Q in Fig. S5.1(c) 6 256 Distance of C to the nearest Pb 4 32 Distance of C to the nearest Pb 4 32 Distance of N to the nearest Pb 4 32 Distance of N to the nearest I 4 32 Table S5.1: 22 features extracted from the internal coordinates. And there are 207 coordinates from 111 system, and 2720 coordinates from 222 system. Some reductant coordinates are removed for the 111 system, like bond angles H1-C1-H2 H1-C1-H3 are kept while H2-C1-H3 are not included for feature bond angle H-C-H, as shown in Fig. S5.1(a). 150 Figure S5.1: (a) MA molecule configuration, shows the bonds C-H, C-N, N-H, angles H-C-N, H- C-H, C-N-H, H-N-H, and dihedral angle H-C-N-H. (b) two octahedra blocks in MAPbI3 system, shows bonds Pb-I, non-bond Pb-Pb, bond angle Pb-I-Pb, I-Pb-I and non-bond angle I-I-I. (c) Two MA molecules shows the relative angles between the orientation of molecules (represented as C-N vector) and axis a and c, center of mass distance between nearby MA molecules. C3. Unsupervised Machine Learning with Mutual Information The Mutual information between variable X and Y are estimated as follow. The support of X and Y is partitioned into bins of finite size based on the K-nearest neighbor (KNN) statistics. 239 Then the mutual information is the summation across these bins as I(X,Y)≈𝐼 binned (X,Y)=∑𝑝 (𝑖,𝑗 ) 𝑖 ,𝑗 log ( 𝑝 (𝑖 ,𝑗 ) 𝑝 𝑥 (𝑖 )𝑝 𝑦 (𝑗 ) ) (S5.7) where i, j are the index for each bin in the support of (X, Y) and p(i,j) = n(i, j)/N is the probability of points falling in this bin (n(i, j), number of points in this bin, N, total number of data points) and px(i) = nx(i)/N, py(j) = ny(j)/N are the probability of points in bin i, j in the 151 support of X and Y, respectively. The entropy of variable X is estimated with the same methodology and is computed as H(𝑋 )≈𝐻 binned (X )=−∑𝑝 𝑥 (𝑖 )log (𝑝 𝑥 (𝑖 )) 𝑖 (S5.8) After the mutual information between the target variable and each individual coordinates or motions is computed, we take the average value from the coordinates for each type of features, which give 44 MI values with NAC (22 from static features as listed in Table S5.1, and 22 from dynamic features), and 22 MI values between static features and the bandgap.

Abstract (if available)

Abstract Nanomaterials exhibit extraordinary thermal, mechanic, optical, electrical and magnetic properties and have wide applications in engineering, materials science and physical chemistry. The theoretical modeling of these nanoscale systems provides a fundamental understanding of various processes, such as physical and chemical synthesis, nuclei and electron dynamics, carrier and energy transfer. Based on the mechanism and timescale, the modeling of these processes requires the development and utilization of theories in classical and quantum dynamics, electronic structure. In this dissertation, two main methods, classical molecular dynamics and non-adiabatic molecular dynamics, will be presented in different scenarios with applied on studying the liquid phase exfoliation process in generating layered materials and the carrier dynamics in carbon nanotube and graphene nanoribbon. We detail the computational schemes along with accurate description for properties of interest. At the same time, we highlight their limitations and disadvantages and the contributions we make to improve their applicability. One contribution we make is to incorporate proper carrier interactions for modeling multi-particle processes with non- adiabatic molecular dynamics. Apart from this, due to the dense manifold of electronic degrees of freedom in nanoscale systems, the modeling requires efficient algorithms and software, and accurate analyzing tools. We incorporate techniques from machine learning into the theoretical modeling on nanoscale systems. We present an unsupervised learning technique for extracting implicit dependency of properties affecting nuclei and electronic dynamics and provides quantitative estimation and comparison of impact from lattice system on the electron recombination. In addition, we present a semi-empirical quantum chemistry model implemented with machine learning framework for improving efficiency and feasibility. The package is designed to take advantage of modern computer architecture graphic process units and is benchmarked with a set of organic molecules. Finally, we present a hybrid model with the developed package and supervised learning technique: artificial neural network. The hybrid model is developed to address the interpretability and transferability problem in machine learning when applying into material science. Incorporating the domain knowledge of quantum chemistry helps alleviate these issues and at the same time, machine learning greatly improves the accuracy of this semi-empirical model. The simplicity, efficiency and accuracy from this hybrid semi-empirical model allows the applications on large molecule systems for modeling long time processes and requiring quantum chemistry level accuracy. This model can be further extended to incorporate with non-adiabatic molecular dynamics for studying processes involving fast and rapid electron dynamics and to push forward the edges of theoretical modeling on nanoscale systems.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Adiabatic and non-adiabatic molecular dynamics in nanoscale systems: theory and applications

PDF

Photoexcitation and nonradiative relaxation in molecular systems: methodology, optical properties, and dynamics

PDF

Dynamics of adiabatic quantum search

PDF

From first principles to machine intelligence: explaining charge carrier behavior in inorganic materials

PDF

Non-adiabatic molecular dynamic simulations of charge carrier dynamics in two-dimentional transition metal dichalcogenides

PDF

Multi-scale quantum dynamics and machine learning for next generation energy applications

PDF

Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials

PDF

Designing data-effective machine learning pipeline in application to physics and material science

PDF

Controlling electronic properties of two-dimensional quantum materials: simulation at the nexus of the classical and quantum computing eras

PDF

Quantum information flow in steganography and the post Markovian master equation

PDF

Towards optimized dynamical error control and algorithms for quantum information processing

PDF

Error correction and quantumness testing of quantum annealing devices

PDF

Trainability, dynamics, and applications of quantum neural networks

PDF

Open-system modeling of quantum annealing: theory and applications

PDF

The effects of divalent counterions on the formation and stabilization of RNA tertiary structure

PDF

Theory and simulation of Hamiltonian open quantum systems

PDF

Theoretical studies of lipid bilayer electroporation using molecular dynamics simulations

PDF

Minor embedding for experimental investigation of a quantum annealer

PDF

Electronically excited and ionized states in condensed phase: theory and applications

PDF

Towards robust dynamical decoupling and high fidelity adiabatic quantum computation

Asset Metadata

Creator Zhou, Guoqing (author)

Core Title Theoretical modeling of nanoscale systems: applications and method development

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Physics

Degree Conferral Date 2021-05

Publication Date 01/13/2022

Defense Date 11/23/2020

Publisher University of Southern California (digital)

Tag Auger process,machine learning,molecular dynamics,molecular simulation,nanomaterials,non-adiabatic molecular dynamics,OAI-PMH Harvest,semi-empirical quantum chemistry

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Boedicker, James< (committee member), Haas, Stephan (committee member), Pilch, Krzysztof (committee member), Wittig, Curt (committee member), Prezhdo, Oleg (dissertation committee chair)

Creator Email guoqingz@usc.edu,kuoching.chou@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-415007

Unique identifier UC11668595

Identifier etd-ZhouGuoqin-9234.pdf (filename),usctheses-c89-415007 (legacy record id)

Legacy Identifier etd-ZhouGuoqin-9234

Dmrecord 415007

Document Type Dissertation

Rights Zhou, Guoqing

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 author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.

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

Repository Email cisadmin@lib.usc.edu