University of Southern California Dissertations and Theses

Neural network for molecular dynamics simulation and design of 2D materials

(USC Thesis Other)

Neural network for molecular dynamics simulation and design of 2D materials

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content NEURAL NETWORK FOR MOLECULAR DYNAMICS SIMULATION AND DESIGN OF 2D MATERIALS By Beibei Wang 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) August 2020 Copyright 2020 Beibei Wang ii TABLE OF CONTENTS Acknowledgements ....................................................................................................................... iii List of Tables .................................................................................................................................. v List of Figures ................................................................................................................................ vi List of Abbreviations ...................................................................................................................... x Abstract ......................................................................................................................................... xi Introduction .................................................................................................................................... 1 Chapter 1: Background ................................................................................................................... 5 1.1 Basics of molecular dynamics ........................................................................... 5 1.2 Dewetting of liquid between MoS2 monolayer ................................................. 8 1.3 Origami 2D metamaterials: graphene ............................................................... 9 1.4 Nanoindentation on Kirigami 2D material: MoS2 .......................................... 11 1.5 Basics of neural networks ............................................................................... 13 1.6 Design of kirigami metamaterials with deep neural network ......................... 15 Chapter 2: Dewetting between MoS2 monolayers ........................................................................ 18 2.1 MD Simulation and validation ........................................................................ 18 2.2 Results and discussion .................................................................................... 20 Chapter 3: Origami structure 2D materials ................................................................................... 27 3.1 Experiment details and MD simulation methods .............................................27 3.2 Results and discussion .....................................................................................29 Chapter 4: Nanoindentation on Kirigami MoS2 ............................................................................ 32 4.1 MD simulation methods ................................................................................. 32 4.2 Results and discussion .................................................................................... 34 Chapter 5: Kirigami material design with deep neural networks .................................................. 40 5.1 MD simulation and machine learning formulation ......................................... 40 5.2 Results and discussion .................................................................................... 45 Chapter 6: Principle of least action for neural network ................................................................ 48 6.1 Action-derived MD with neural network ........................................................ 48 6.2 Neural network structure and simulation methods .......................................... 51 6.3 Results snd discussion .................................................................................... 54 Chapter 7: Summary and future directions ................................................................................... 59 7.1 Summary ......................................................................................................... 59 7.2 Future work ..................................................................................................... 61 References .................................................................................................................................... 63 iii Acknowledgements First of all, I would like to express my sincerest gratitude to my Ph.D. advisor, Prof. Rajiv Kalia, for his continuous care, support and valuable suggestions over the past many years. Under his guidance, I not only learned how to perform molecular dynamics simulations, get insights from data and prepare high-quality academic manuscripts, I also had many great opportunities to widen my perspective for innovative research and to enhance my skills for challenging unsolved scientific problems. Prof. Kalia has taught me how to become a successful researcher in every aspect, his enthusiasm for research, patience for students and his immense knowledge has inspired me throughout my Ph.D. career. I could not imagine having a better Ph.D. advisor and I will forever be thankful to Prof. Kalia. Besides my advisor, I express deep gratitude to Prof. Priya Vashishta and Prof. Aiichiro Nakano. I benefited greatly from their enlightening courses in classical and quantum molecular dynamics. Their generous support and care, their earnest attitude in research and rigorousness in work are inspiring. And I will always remember Priya’s lessons, such as “complex begins simple”, “never make a decision when you are emotional”, “it is fine to admit a lack of knowledge” and so on. I will also remember the 2017 Chinese New Year with Aiichiro and his family. Priya and Aiichiro are my best role models for a teacher, researcher and mentor. I’m grateful to other faculty members on my Ph.D. committee. I thank Prof. Satish Kumar for insightful discussions on related research and I thank Prof. Stephan Haas for suggestions and help in general. Their presence as committee members made this defense possible. I would also like to thank my research collaborators. Dr. Pankaj Rajak is Prof. Vashishta’s former student, I thank him for fruitful discussions over many research topics and I hope that I could be as enthusiastic and energetic as him. Dr. Huan Zhao, who is a former student at USC, I thank him for collaborating with me on research about origami 2D materials. I thank Dr. Mark iv Stevens, a scientist from Sandia National Lab, for suggestions and help in research and contributing to some of my manuscripts. I would like to acknowledge the support from the U.S. Department of Energy (DOE). My research projects were funded by the DOE Office of Science, Basic Energy Sciences, Materials Science and Engineering Division. Some of my work was performed at the Center for Integrated Nanotechnologies, an Office of Science User Facility operated for the DOE Office of Science. Last but not least, I would like to thank my family and friends. I thank my father, my mother and my sister for trusting and supporting me throughout writing this thesis; I thank my parents in law and my wife for taking good care of my son and in general; I also thank all my friends (too many to list here) for offering support and friendship that I needed over the years and for the many years to come. v List of Tables Table 2.1: Parameters for Coulomb and L-J interactions between H2O and MoS2 P-23 Table 4.1: Parameters for L-J interactions between MoS2 and the indenters P-36 Table 4.2: Mechanical properties for H patterns P-41 Table 4.3: Mechanical properties for S patterns P-41 Table 5.1: DQN training pseudocode with experienced replay P-45 Table 6.1: L-J interaction parameters and lattice constant P-52 vi List of Figures Figure 1.1: a simplified illustration of a standard MD simulation process. P-14 Figure 1.2: A schematic of a neural network. The input layer (black) consists of two units, taking in x1 and x2, respectively. The hidden layer (blue) has 4 neurons and the output layer (red) gives an output vector (y1, y2, y3). The blue and red edges connecting layers represent the weights, each node also applies a bias to its input. P-20 Figure 2.1: MD simulation schematic set up for contact angle simulation. O and H atoms are red and white, respectively; S atoms are yellow, and Mo atoms are not shown in here. P-24 Figure 2.2: Density profile of a H2O droplet and the contact angle as a function of droplet size. (a) schematic diagram of the contact-angle measurement. The colored map represents the density distribution in cylindrical coordinates. Z is the height from MoS2 substrate and r is the distance from the central axis of the droplet. ! is the contact angle at the interface of H2O and MoS2. (b) contact angle as a function of radius of the contact area. P-24 Figure 2.3: MD simulation set ups for (a) H2O monolayer, and (b) IPA/H2O mixture between MoS2 bilayers. For visualization purposes, MoS2 layers are lifted to expose the liquid layer. Oxygen and hydrogen atoms are red and white, Mo and S atoms are pink and yellow, respectively. IPA molecules are brown. Periodic boundary conditions are applied in the x and y directions of MoS2 bilayers. The size of the system is ~100 ´ 100nm 2 in the x-y plane. P-25 Figure 2.4: Dewetting of an H2O film between MoS2 membranes. Red regions represent H2O. (a) Top view of the H2O monolayer at t = 0 ps. Dry patches appear spontaneously and start expanding rapidly. (b) H2O molecules form a network structure at t = 50 ps. (c) At t = 200 ps, the network breaks up into H2O droplets. Some of them coalesce to form larger droplets. (d) H2O dewetting at t = 500 ps. (e) Radius of a dry hole in the film expands linearly with time and the slope gives a dewetting speed of 373 m/s. P-26 Figure 2.5: Temperature and pressure of the system during dewetting. (a) “Temperature” of water rises quickly in the first 50 ps of dewetting. Subsequently, the energy transferred from water to MoS2 increases the temperature of MoS2 membranes. The equilibrium is reached after 300 ps. (b) Shows that pressure increases rapidly as the water monolayer breaks up into nanodroplets. P-27 Figure 2.6: Dewetting of an IPA/H2O mixture between MoS2 membranes. (a) Snapshot of the IPA/H2O liquid film t = 0 ps. (b) Snapshot of the mixture taken at t = 0.5 ns shows a percolation network in which IPA and H2O molecules have phase separated. IPA molecules are mostly outside and H2O are inside the network. (c) Snapshot of the IPA/H2O mixture t = 2.0 ns shows an increase in the fraction of dry patches. (d) shows the formation of interconnected nanodroplets after 10 ns. (e) Shows a linear increase in the radius of a circular dry patch (inset). The speed of the dry patch in the mixture is significantly less than that in pure water. P-28 Figure 2.7: Temperature and pressure of the IPA/H2O film during dewetting. The temperature (a) and pressure (b) of IPO/H2O change rapidly in the first 500 ps. Subsequently, the temperature of the mixture and MoS2 membranes change slowly and thermal equilibrium is not established even after 3 ns. In contrast, the equilibrium is reached within 200 ps in the absence of IPA. P-28 Figure 2.8: Time evolution of dry patches in H2O and 50% IPA/H2O mixture during dewetting between MoS2 membranes. Dewetting is much slower in the mixture case because of the slow diffusivity of IPA molecules. The fraction of dry patch in the mixture is lower than that of pure H2O, indicating that the mixture covers more surface area of MoS2. P-29 vii Figure 2.9: Deformation of MoS2 membranes caused by dewetting of the H2O film. (a) Dimples are formed in MoS2 membranes after dewetting. (b) Close-up view of deformation in MoS2. (c) Shows a strip cut from the deformed MoS2 cap. Knowing the bending energy and bending force on the strip, we calculate the bending modulus to be 100±20 N/m. P-29 Figure 2.10: Water molecules trapped between MoS2 membranes form two types of structures. (MoS2 membranes are lifted to visualize the structures of water.) (a) Water nanodroplet confined inside an MoS2 cage. (b) Oxygen-Oxygen radial distribution function of H2O in the confined nanodroplet. (c) Confined water molecules form a monolayer of “frozen” H2O in registry with the MoS2 lattice, and (d) shows the triangular lattice structure formed by water molecules. P-30 Figure 3.1: Demonstration of the folding technique. a) A step-by-step illustration of the folding technique; b) A few-layer graphene flake on PVA/SiO2/Si substrate before and after folding. The AFM image inset shows the folding crease is well aligned with the user-defined PVA edge. c) Deterministic folding of patterned CVD graphene monolayer array under different scales. Dashed and solid arrows indicate the direction of the triangle tip before and after folding, respectively. There is a rectangular PVA layer under each triangular graphene flake, which can be observed by trained eyes. d) An AFM image of a folded graphene sample (after PVA layer is removed) showing the high- quality interface, as well as the well-defined folding edge and folding angle. The inset is the optical microscope image of the same sample, where the area for taking AFM imaging is outlined by a dash line. P-31 Figure 3.2: MD simulation setup. (a) The simulation box is full of water molecules and in the middle there is a triangular graphene sheet. White spheres are hydrogen, red is oxygen and cyan is carbon atoms. Only half of water is shown in the figure for visualization purpose. We apply constant pressure P = 0.08 GPa on water from right side of the box as the yellow arrow indicates. (b) Top view of graphene sheet in the system. Dashed lines are the boundaries between harmonically-fixed and free-to-move C atoms. The red line is parallel to the left edge of the graphene sheet, the black line is a “tilted” boundary. # is the orientation angle of the boundary with respect to the left edge. In the simulation we have chosen # = 0˚ and 12˚. P-32 Figure 3.3: folding angle of graphene as a function of simulation time. (a) The graphene patch is folded by fluidic force gradually. (b), (c) and (d) are side views of the simulated water (red and yellow) and graphene (cyan) system with a graphene bending angle of 90º, 120º and 180º, respectively. Note, the region colored with yellow is the fraction of water molecules to which we apply a constant pressure to drive the fluidic flow. P-33 Figure 3.4: Understanding the folding mechanism by molecular dynamics simulation. Snapshots (a), (b), (c), and (d) show graphene bent at angles of 0º, 90º, 150º, and 180º, respectively. Here, oxygen and hydrogen atoms in H2O and carbon atoms in graphene are red, white and cyan spheres, respectively. Yellow spheres represent molecules in the water wake. Note that only half of the water in the MD box is shown in order to visualize the folding of graphene. (e) shows accumulated work per atom done on graphene by water flow as a function of time. Note that φ represents the orientation angle of the boundary between fixed and free portions of graphene with respect to the y axis of the simulation box. (f) shows changes in the internal energy of graphene as a function of the folding angle. The internal energy increases, reaching a maximum value for q = 120º. The potential energy barrier for graphene folding is around 2 to 3 meV. P-34 Figure 4.1: A schematic of nanoindentation simulation (a) Top view of the centrosymmetric nested hexagonal kirigami pattern in a MoS2 monolayer. (b) Top view of square pattern in MoS2. In both (a) and (b), the white region are the places where we have removed molybdenum and sulfur atoms. Yellow region represents the sulfur atoms. (c) Snapshot of the conical diamond indenter we used in the nanoindentation simulation. Cyan spheres are carbon atoms. The indenter angle at the tip is 120 degrees. (d) and (e) are snapshots of conical indenter on top of the MoS2 membrane. P-35 viii Figure 4.2: (a) Load-displacement curves for nanoindentation on MoS2 with 3 kirigami structures of different area fractions. As a reference, the response of the MoS2 monolayer without the kirgami structure is shown in fig. S2 in the supplementary material. The green curve is the indentation response curve of the kirigami pattern shown in (b). The blue and red curves correspond to patterns shown in (c) and (d). The arrows in (a) pointing upward and downward indicate loading and unloading processes, respectively. The hysteresis in load-displacement curves indicates irreversible damage in MoS2 kirigami structures during nanoindentaion. The MoS2 kirigami pattern in (d) can be indented much more with a smaller load than the pristine MoS2 monolayer. P-37 Figure 4.3: (a) Load-displacement curves for nanoindentation on square MoS2 kirigami patterns shown (b), (c) and (d). The upward and downward arrows in each curve indicate loading and unloading processes, respectively. The area of kirigmai patterns in (b), (c), and (d) are 5.4%, 13.3%, and 24.7%, respectively, of the total area of the MoS2 monolayer. As in the case of hexagonal kirigami patterns, the hysteresis and hardness decrease significantly with an increase in the kirigami area. P-38 Figure 4.4: Effect of the indenter size on the two MoS2 kirigami structures. (a) Load vs indentation depth for the hexagonal kirigami structure indented with indenters of radii 30 Å, 45 Å and 60 Å. (b) Load vs indentation depth curves for the square kirigami structure indented with the same set of indenters. The area fractions of kirigami patterns in panels (a) and (b) are ~ 25%. P-39 Figure 4.5: Stress distribution in the deformed MoS2 membrane. (a) and (c) are side views of the deformed MoS2 kirigami structures. MoS2 layers form an elastic pyramid-spring under nanoindentation. (Indenter is removed from the view for better visualization of MoS2). Stress is accumulated symmetrically around the corners and arms of kirigami patterns and goes up to 10 GPa. Panels (b) and (d) show top views of deformed MoS2 under nanoindentation. (Stress is calculated using per-atom volume and per-atom stress tensor averaged over 0.2 million MD steps.) P-40 Figure 4.6: Evolution of defects induced by nanoindentation in MoS2 kirigami structures. (a) close- up views of a corner in a hexagonal kirigami pattern before indentation (left panel) and after indentation (middle panel). There are broken bonds and a nanocrack at the corner after the indenter is removed and the system is relaxed. (b) shows irreversible damage caused by indentation in a square kirgami structure. The left panel shows the atomic configuration at a corner of the square kirgami pattern before indentation. The middle panel is a close-up view of the same corner after indentation. Point defects shown in the middle panel persist even after the indenter has completely retracted and the kirigami structure has been relaxed for 2 million MD steps. The defect density in the two systems is shown in (c). Blue curves represent defect density in the hexagonal kirigami structure in figure 1(d) during loading (lower curve) and unloading (flat upper curve) phases. Red curves in (c) show the defect density during loading and unloading in the square kirigami structure in figure 4.3(d). P-41 Figure 5.1: illustration of MoS2 kirigami structures and their corresponding stress-strain curves. (a) the horizontal lines divide the MoS2 nanosheet into independent 6 rows, in each row we define 13 different configurations of a cut slit. All cuts in the 6 rows combined form a kirigami pattern in the MoS2 nanosheet. (b) The stress-strain curve corresponding to the kirigami structure in (a). (c) examples of possible kirigami patterns in our design space. P-43 Figure 5.2: An illustration of the RL model. (a) a schematic of the RL agent, the input to the agent is a binary tensor representing the kirigami structures and the output is the corresponding stretchability, i.e. the maximum strain that this kirigami MoS2 can withstand. (b) loss function against number of training epochs. (c) the expected reward as a function of training epochs. P-44 ix Figure 5.3: distribution of max strain and comparison between RL-generated and randomly selected structures. (a) shows the distribution of the maximum strains for the entire state space for n = 4 and sampled space for n = 6. For both cases, the majority of the structures have max strain less than 20%. (b) distribution of the max strain of 500 RL generated kirigami structures, along with a visualization of some of the structures. (c) distribution of 500 kirigami structures uniformly sampled from the database, and some typical structures with their strechability. P-46 Figure 5.4: (a) Principle component of the state-action value function against cut sequence. (b) illustration of the creation of kirigmai structure by the RL agent. P-47 Figure 5.5: atomic stress distribution and the stress-strain cuve. (a) sideview of a brittle MoS2 kirigami structure that is about to fracture, no out-of-plane deformation is observed. (b) and (c) different components of the atomic stress, tensile stress accumulates at the weak points of the structure. (d) side view of a ductile MoS2 structure, where out-of-plane deformation happened. (e) and (f) not only tensile stress is concentrating at the corners and bridges of the cut pattern, we also observe compressive stress where out-of-plane bending is observed. (g) the stress-strain curve of the two structures shown in (a) and (d). P-48 Figure 6.1: A schematic of the neural network. The system configuration for all spatial degrees of freedom is calculated by inserting an instant of time t into the input layer. A single hidden layer (blue) is used and the output layer (red) gives the atomic coordinates whose first and second derivatives provide velocity and acceleration, respectively. The blue and red lines connecting layers represent the weights and each node also has a bias. P-51 Figure 6.2: Four selected NN (solid blue) and MD (solid red) trajectories show that the NN can reproduce non-trivial paths. These are for a time period of 1 ps. P-54 Figure 6.3: Comparison between the neural network (NN) and MD simulation. (a) Kinetic (blue), potential (green) and total (red) energy per particle as a function of time for MD (square dots) and NN (solid lines). Radial distribution function and velocity autocorrelation function for MD (blue) and NN (red) are shown in (b) and (c), respectively. P-55 Figure 6.4: (a) shows how the three components of the loss function decrease with the number of epochs. The boundary condition term (red) and energy conservation term (blue) fall more rapidly than the OM action term (green). (b) Shows the mean square error (MSE) for trajectory (red) and momentum (blue) which was calculated by comparing the MD trajectory to the neural network predicted trajectory. (c) shows differences between MD and NN trajectories (red) and velocities (blue) as a function of time. Here t = 2ps. P-56 Figure 6.5: Scaling of NN for 108-particle system. (a), (b) and (c) are the run time for every 1000 training epoch as a function of the number of units in the hidden layer, the size of the time grid, and the number of atoms in the system, respectively. P-56 x List of abbreviations LM layered material TMD transition metal dichalcogenides LPE liquid-phase exfoliation MD molecular dynamics ANN artificial neural network PDE partial differential equation LJ Lennard-Jones tBLG twisted bilayer graphene AFM atomic force microscopy EBL electron beam lithography REBO reactive empirical bond order IPA isopropanol alchohol PVA polyvinyl alcohol CNN convolutional neural network DQN deep Q-network xi Abstract In recent years, two-dimensional (2D) materials have become the spotlight of nanoscience owing to their exotic properties. Massive experimental and theoretical research has been triggered to understand the synthesis of 2D materials via exfoliation and to explore the controllability of their properties via structure engineering. As a synthesis technique, liquid-phase exfoliation is of great significance due to its simplicity and high yield. Researchers have devoted many efforts to explain and improve the synthesis of 2D materials. However, understanding and optimizing the synthesis experiment at nanoscale remains challenging. Structure engineering methods such as origami and kirigami offer efficient approaches to designing and tuning novel 2D metamaterial structures and properties. Folding and stacking of atomically-thin nanosheets can lead to intriguing new physical phenomena such as bandgap tuning, metal-insulator transition, Van Hove singularity and superconductivity. Nevertheless, achieving scalable folding and reasonable designing of 2D materials with high spatial precision and desired properties have been difficult tasks. To better understand the synthesis of 2D materials via liquid-phase exfoliation, to develop a general method for origami and kirigam structure designing, we first perform molecular dynamics simulations to study the interlayer coupling and delamination of a 2D material in the presence of liquid solvents and then propose a scalable deterministic technique to fold the synthesized 2D metamaterials with a pre-defined configuration. We investigate the mechanical properties of kirigami structures of a 2D material (i.e. MoS2) and discuss the underlying mechanism and energetics of the process. Finally, we demonstrate how neural-network based models can accelerate the design of kirigami 2D materials and address foundational problems in molecular dynamics. 1 Introduction Since the experimental synthesis of monolayer graphene by mechanical exfoliation of bulk graphite[1], extensive research has been devoted to the field of two-dimensional (2D) layered materials(LMs). Many new 2D materials have been discovered or re-discovered, such as transition metal dichalcogenides (TMDs)[2], boron-nitride (h-BN)[3], phosphorene[4], [5] and so on. The family of 2D materials now offers a wide range of electronic properties, from conducting graphene to semiconducting molybdenum disulfide (MoS2) and insulating h-BN. Moreover, the 2D crystal structures leverage a unique combination of mechanical properties, with high in-plane stiffness and strength but extremely low out-of-plane flexibility[6]–[8]. Paired with these novel properties, 2D materials are promising candidates for a wide range of applications in various fields such as catalysis[9], photonics[10] and spintronic devices[11]. As an example, researchers have created sensors and high-performance electrodes using MoS2[12]. Their outstanding properties and the controllability of these characteristics are becoming the spotlight of research in Nanoscience. Among all the 2D materials, transition metal dichalcogenides (TMDs) are drawing significant attention in research and is considered as a worthy candidate for next-generation thermoelectric and optoelectronic devices due to their excellent electrical properties[13], mechanical strength[14] and chemical stability[15]. For example, MoS2 is known to be a solid lubricant and a promising new electronic material for wearable electronic devices due to its high electron mobility, wide band gap[16] and excellent flexibility[17]. In fact, electronic and mechanical properties of these TMD’s can be easily regulated with artificial defects, mechanical straining or by creating hybrid material consisting of vertically stacked heterogenous TMDs. For example, Apte et al. have shown experimentally that under strain MoWSe2 heterostructure goes structural transformation near crack tip from semi-conducting 2H phase to 1T metallic phase, which increases its fracture toughness[18]. 2 To synthesize TMDs and other 2D materials, one common mechanical approach is liquid- phase exfoliation (LPE) by sonication [19]. Take MoS2 as an example, a recent MD simulation reveals that, in sonication experiment, tiny bubbles collapse at the interface between bulk MoS2 and the solvent, creating nanojets of liquid molecules that strike the bulk material with a high velocity [20]. The interaction between the bulk and the shockwave from nanojets enlarges the gap between MoS2 layers, then liquid molecules get between the layers, weakening the intra-layer coupling. However, a complete dewetting process can expel these intercalated solvent molecules due to the re-combination of MoS2 layers by intra-layer coupling. Therefore, understanding dewetting behavior of the liquid layer between MoS2 interfaces is important for the synthesis of MoS2 and other TMDs by liquid-phase exfoliation. In this thesis, we will report dewetting of different solvents between MoS2 layers using MD simulation. With the development of efficient experimental preparation of layered materials, many innovative designs of metamaterials have been proposed experimentally and theoretically. For example, Hwang et al. demonstrated stretchable metamaterial electronics and kirigami actuators with highly tunable mechanical responses by introducing major and minor cuts to control kirigami deformation[21]. Tang et al. designed programmable kirigami-based metamaterials and showed that the tilt orientations of cut units can be controlled. They also simulated the application of this design in reducing building energy consumption[22]. Despite the importance of kirigami, it is a challenge to rationally design artificially structured metamaterials with desired mechanical and physical properties that are not accessible in nature materials. In this thesis, we report MD simulations to study two typical metamaterial structures: origami and kirigami. As an ancient paper art, origami has recently become an emerging technique for fabricating LM based devices. Graphene as a 2D material, exhibits exceptional mechanical properties such as ultra-flexibility, high strength and low bending stiffness[6], [7], [23]. Owing to these extraordinary 3 properties, graphene can be folded into various structures [24], leading to exceptional physical properties leveraged by the strong intralayer interaction. Recent work about twisted bilayer graphene (tBLG) has demonstrated fundamentally intriguing electronic properties such as Mott insulator transition[25], Van Hove singularity and superconductivity[26]–[28], which shed light towards Moiré physics and twistronics[29]. As these properties are dictated by precisely tuning the interlayer misorientations and misalignment, a simple technique that can create twisted 2D materials with high controllability and scalability is desirable. In this thesis, we will describe an experimental method for folding 2D materials with a pre-defined orientation using micro-fluidic forces. This technique is easily controllable and can twist 2D materials without the need of manually align and stack the material layers. Experimental details can be found in the work by Zhao et al.[30] To understand the mechanism and energetics of the folding, we perform MD simulations of folding of graphene under fluidic flow, which will be described in chapter 2.3. As mentioned above, another way of creating complex 3D structures from 2D material is kirigami. Unlike origami, kirigami is an art of paper-cutting. 2D layered materials such as graphene and TMDs are known to be as flexible as paper. To explore their flexibility, researchers have recently applied kirigami technique to these 2D materials and made 3D metamaterials with exceptional mechanical and electrical properties. For example, Blees et al. created graphene-based devices with two types of kirigami patterns: pyramid spiral spring and alternating offset cuts[7]. Their measurements of mechanical and electrical properties indicate the flexibility and viability of graphene-based flexible electronic devices. In this thesis, we will study the flexibility and stretchability of MoS2 kirigami structures and guide the creation of MoS2 metamaterials with kirigami pattern to maximize the yield strain under constraints. In addition to MD simulations, we also combine machine learning (ML) method with MD. We have studied the mechanical properties 2D MoS2 metamaterials using reinforcement learning 4 (RL) and applied an artificial neural network (ANN) approach in foundational MD research. In recent years, ANN has become very popular not only in computer science, but also in nanoscience. Shimamura et al. developed an ANN-based empirical interatomic force field using data from first- principles MD simulations[31]. Fukushima et al obtained the free energy of alkali metals at low computational cost using an ANN force field[32]. NN has also been applied to solve differential equations. Piscopo et al demonstrated a flexible approach to solve the second-order partial differential equations (PDEs) in cosmology and studied phase transitions of the early universe[33]. In this thesis, we describe two scenarios where ML approach is applicable: PDE solver based on the principle of least action, and RL-assisted design of kirigami metamaterials. The thesis is organized as follows: in chapter 2 we will introduce the background of MD simulation, origami and kirigami structure of 2D metamaterials, then the application of neural network in MD simulation and design of 2D materials. In chapter 3, more details about our simulation study on the intra-layer coupling of MoS2 is discussed. Chapter 4 and 5 will focus on the origami and kirigami design of 2D MoS2, respectively. In chapter 6, we use a simple case to illustrate the application of ANN and the principle of least action to solve a boundary-value problem in classical MD simulation. Chapter 7 will further discuss the application of ML in designing 2D metamaterials. Chapter 8 is a summary of the current work and an outlook for future directions. 5 Chapter 1: Background 1.1 Basics of molecular dynamics MD is a simulation method for analyzing the physical state of atoms and molecules at the microscale. Given the form of interatomic potential, we use MD to mimic the movement of atoms based on Newton’s second law: ) ⃗ = +,⃗. In a many-particle system, atoms interact with each other for a certain time period, usually nanoseconds or microseconds, depending on the time scale of the interested physics process. MD allows us to have a view of the dynamic development of the system and get insight on how underline principle of physics governs the evolution of and interaction between different components of the system. The basic idea behind MD is to find the solution to Newton’s equation of motion. In non- trivial cases, the movement of all the particles are coupled, therefore it is a many-body problem and the algebraic solution is not feasible. However, a numerical solution is mathematically straightforward using finite difference schema: - ⃗(/+∆/) = - ⃗(/)+3 ⃗(/)∆/ (1.1) 3 ⃗(/+∆/)= 3 ⃗(/)+,⃗(/)∆/ (1.2) where - ⃗(/), 3 ⃗(/) and ,⃗(/) are position, velocity and acceleration of a particle at moment t, respectively. The idea is to integrate the equation of motion simultaneously for all particles in the system. In practice, we often use the so-called “time symmetric” and more stable integration methods such as “velocity-Verlet” algorithm[34]: - ⃗(/+∆/) = - ⃗(/)+3 ⃗(/)∆/+ 1 2 ,⃗(/)∆/ 4 (1.3) 3 ⃗(/+∆/)= 3 ⃗(/)+ 1 2 [,⃗(/)+,⃗(/+∆/)]∆/ (1.4) 6 where acceleration ,⃗(/) and ,⃗(/+∆/) of particles are calculated based on the coordinates - ⃗(/) and - ⃗(/+∆/) of every particle in the system using ) ⃗ = −∆8({- ⃗}) and Newtons’ second law. A schematic illustration of a typical MD algorithm can be found in the following figure 1.1. Figure 1.1: a simplified illustration of a standard MD simulation process. MD is a powerful tool. Many research works have been done to explore physics process in the world at nanoscale and microscale using MD. With this method, researchers are able to test theories of interatomic and molecular interactions and develop interaction models to describe chemical reactions[35]; MD can also help us interpret the results of experiments[30], identify novel structures and dynamics of nanomaterials and therefore propose applications of the materials. MD improves our understanding of fundamental physics through modeling[36]. However, MD is not the elixir, there are limitations of classical MD. As an example, one of them is the accuracy issue. Compared with quantum description of atomic and molecular interactions, MD force fields are inherently approximations. For instance, a typical force field often incorporates Lennard-Jones (L-J) terms that capture short-range repulsive and long-range attractive interactions between atoms; spring-like interactions model the equilibrium length of covalent bonds, etc. To tackle different research problems, one needs some experience to choose 7 a proper force field. Ho et al performed MD simulations to study the graphene-water interface with multiple force fields for water[37]. They found the TIP5P[38] model deviates the most from other models in terms of the interface structure of water, and they observed a large variation in the predictions of dynamical properties for water. In this thesis, we have made the first-step effort to address two foundational limitations of conventional MD method. One is the timescale issue. MD simulations are usually performed in the domain of nanoseconds and rarely microseconds[39], where many physics processes happen, e.g., folding of proteins[40], [41]. The other issue with MD is the search space complexity. One important research topic in metamaterials is the design of materials themselves. Screening the structures and compositions based on material properties via brute force ab initio simulation is not only expensive but also exceeding the limits of the modern computing resources. Many efforts have been devoted to address these two limitations. Nakano designed a scalable parallel algorithm based on the elastic band method and studied long-time dynamics of many-particle system[42]. Bassman et al. developed a platform for accelerating the design of layered materials using Gaussian process regression model and Bayesian optimization method[43]. In this thesis, we will introduce two ML methods as our attempts to realize time parallelization of classical MD simulation and assist the design of kirigami metamaterials, respectively. There are several packages to run MD simulations, just to name a few: GROMACS (GROningen MAchine for Chemical Simulations), Amber, NAMD (NAnoscale Molecular Dynamics), Desmond, OpenMM, etc. A very popular package is called LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator)[44]. It is a general-purpose simulation package that runs very fast on parallel-computing machines. It has detailed instructions and examples, easy to modify and extend. LAMMPS read customized control commands from a input file, where we initialize system configuration, set up force field parameters, choose statistical ensemble and 8 proper timestep. There are customizable commands that allow users to control temperature, pressure the deformation of the system if needed. All the MD simulations introduced in this dissertation is performed using LAMMPS. 1.2 Dewetting of liquid between MoS2 monolayer MoS2 as one of the spotlight 2D materials has been extensively reported by researchers using MD [45]–[48]. The interlayer coupling between Mo-Mo, Mo-S and S-S atom pairs can be described using van der Waals interaction[49]. However, only a limited number of simulations have been performed to investigate the interaction between MoS2 and adsorbates[50]. The motion of nanofluidic molecules at MoS2 interface is of great importance, it leverages the potential applications of MoS2 nanostructure in biological sensing[51], molecular filtration and DNA sequencing[52]. Structure and dynamics of liquid molecules at the nanoscale is under extensive investigation, many simulations have been performed by confining water molecules and clusters in nanocontainers[53], nanoporous structures[54] and so on. Gao et al. studied alcohol/water mixture confined under graphene slit with MD[55]. They evaluated the four types of mixtures under confinement to evaluate the effect of Alcohol layer on the dynamics of water molecules and proposed potential application in exfoliation of graphene. Understanding the behavior of liquid films at solid interface is important to a variety of applications such as self-removal of liquids, antifogging[56], self-assembly of microscopic clusters in liquid[57] and water harvesting technology[58] in which liquid molecules on the solid surface aggregate and condensed into droplets through dewetting. Many experimental studies and MD simulations have been carried out to investigate the dewetting behavior of liquids confined to solid interface. Jensen et al. have studied the role of interfacial energy in dewetting of water on hydrophobic surface[59]. Kayal et al. used MD to study dewetting of water in carbon nanotubes, they found the tunnel flow and orientation of water molecules can be tuned by an electric field 9 applied normal to the direction of transport[60]. Zhang et al.’s MD study of water nanofilms in contact with silica surface indicates that dewetting is a two-stage process in which dry patches appear due to thermal fluctuations and water film contracts because of hydrogen bonding and electrostatic interactions between water and the surface[61]. Dewetting is an important phenomenon in the synthesis of 2D materials. In MoS2 LPE experiment, solvent molecules enter into the intra-layer space of the bulk material under the shear stress on the material, weakening the inter-layer coupling of MoS2[62]. One commonly used organic solvent for this experiment is isopropanol and water mixture (IPA/H2O)[19], [63]–[65]. In this thesis, we report MD simulation studies to examine the dewetting behavior of pure H2O and H2O/IPA mixture confined between molecularly thin MoS2 membranes. Our results reveal distinct dewetting processes of water and the mixture in the galleries of MoS2 bilayers and provide explanation for the effectiveness of the mixture over pure water in LPE experiments. 1.3 Origami 2D metamaterials: graphene Graphene and other atomically thin 2D materials have demonstrated paper-like mechanical characters such as low bending rigidity, high flexibility and strength[6], [7], [23]. With this exceptional properties, graphene can be folded into origami structures, giving rise to exotic physical properties and various applications[66], [67]. Through systematic MD simulations, Zhu et al. proposed a robust hydrogenation-assisted method to realize a programmable and reversible fabrication process for graphene origami nanocage, and demonstrated a promising application in molecular mass capture, storage and release[24]. Inspired by origami, Mu et al. fabricated a graphene-paper based walking device that can adopt predesigned shapes and motions through self- folding mechanism[68]. Recently, Cao et al investigated the electronic behaviors of magic-angle bilayer graphene superlattices[25]. They fabricated the tBLG structures by stacking two sheets of graphene and twist 10 them relative to each other by a small angle, then observed unconventional superconductivity with critical temperature of up to 1.7 Kelvin[28]. However, consecutively transferring two graphene monolayers and precisely aligning them with the magic angle is difficult, since it is hard to identify the crystal orientation without sophisticated characterization tools[69], therefore it offers limited controllability over the angle and moderate scalability with this technique. In this thesis we report a universal experimental method for twisting and folding 2D materials with a predefined position and orientation using micro-fluidic forces. This technique can easily prepare highly controllable tBLG and other twisted 2D materials without the need of manually aligning and stacking the material flakes, nor characterizing the crystal orientation. Using this method, both micrometer-scale exfoliated material flakes and wafer-scale CVD-grown 2D materials can be folded in seconds, with sub-micrometer spatial accuracy. Starting from 3% wt. polyvinyl alcohol (PVA) water solution, we spin-coated it onto Si substrate, followed by a two- minute backing at 90℃ to form a uniform PVA layer with a thickness of ~40 nm confirmed by atomic force microscopy (AFM), then this 2D film was transferred onto PVA/Si substrate with mechanical exfoliation method. The next step is to expose the designed sample area by a standard electron beam lithography (EBL) technique, to cross-link the PVA with a dose of ~ 5000 µC/cm 2 . After EBL, the sample was soaked into 50℃ deionized (DI) water. A water flow is applied over the substrate with a speed of ~10 µm/s, and the direction of water flow is normal to the pattern crease to assist the folding. Experimental details can be found in the work[30] by Zhang et al. To understand the physical mechanism underlying the folding process, we carried out MD simulation to model the folding of monolayer graphene. We monitored the folding angles of graphene, investigated the free energy change and examined the effect of hydrophobic graphene on the water flow for two different folding configurations of graphene nanosheet. A detailed description of simulation setup and results is given in chapter 3. 11 1.4 Nanoindentation on Kirigami 2D material: MoS2 MoS2 as a 2D layered material, although recently have attracted great attention in research for flexible electronic devices owning to its unique structure, novel electronic properties[70] and high mechanical strength[14], tensile stretchability of MoS2 and other 2D materials limits their application in realizing fully foldable, bendable and stretchable devices[71][72]. Efforts have been made to integrate super elasticity in MoS2 by introducing kirigami. Unlike paper folding, kirigami is an art of paper cutting. Known to be as flexible as paper, 2D materials such as graphene and TMDs have been explored by researchers using kirigami technique[7], [73]. Structures with exceptional mechanical and electrical properties have been fabricated. For example, Zheng et al. designed electronic architectures of MoS2 at nanoscale[74], and they found adopting kirigami patterns can greatly improve the stretchability and electrical properties for 2D MoS2-based devices. The flexibility of 2D materials is best characterized by the Föppl-von Kármán number: ; = < => ? = @ , where A 4B represents the 2D Young’s modulus of the material, C is the out-of-plane bending stiffness and D 4 is the area of a square sheet of the material[75], [76]. Blees et al. performed experiments to obtain ; for graphene, and they found a large value of ; in the range of 10 5 -10 7 [7]. We have estimated ; for a MoS2 monolayer using the value of C and A 4B calculated from MD simulation with the active empirical bond-order (REBO) force field[49], [77]. With these values, the ; number for a 100´100 nm 2 sheet of MoS2 is found to be ~1.2´10 6 , which falls in the region of ; for graphene and close to that of a standard paper sheet. Such a large value of ; indicates that kirigami patterns can be designed experimentally in a MoS2 monolayer membrane. As an important experimental technique to study the hardness and mechanical properties of materials at the nanoscale, nanoindentation has been applied to 2D materials widely[23], [78]– [80]. Bertolazzi et al. performed nanoindetation experiment to measure the mechanical strength of 12 single and bilayer MoS2[14]. They have found that the in-plane stiffness and breaking strength of ultrathin MoS2 are comparable to those of stainless steel. Stewart et al. carried out MD simulations to study nanoindentation on the basal plane of thin MoS2 layers[81]. They studied the incipient plastic deformation and found the local phase transformation using a slip vector analysis. Tan et al. performed simulation to study the deflection dependency of monolayer graphene and obtained modulus from indentation[82]. Recently, some simulations and experimental studies of 2D kirigami materials have been reported. Shyu et al. reported a controllable kirigami patterning technique to engineer the elasticity of nanocomposite sheets[83]. They claimed to have significantly improved the tensile flexibility and the electrical properties of the kirigami-based stretchable conductor graphene nanocomposites. Hanakata et al. performed simulations to investigate the density and length of overlapping kirigami patterns in MoS2 nanoribbon[84]. They found that kirigami technique significantly enhanced the tensile ductility and the yield strain of MoS2. In chapter 4, we report systematic MD simulations to study the elastic behavior of two differently patterned kirigami MoS2 monolayers under nanoindentation. We use the previously mentioned REBO potential to model the interactions for our Mo-S system. The REBO potential accounts for changes in local atomic configurations of atoms which allows us to analyze broken bonds and defects. The binding energy of our system can be expressed as: E F = 1 2 G) HI J K- HI L[M1+ N - HI OPQ RST UV −W HI XQ RYT UV ] HZI (1.5) where ) HI J is a cut-off function, - HI is the distance between atom pair [ and \, W HI is a bond-order function, and P, X, N, ] and ^ are parameters for pairwise interactions[49]. 13 1.5 Basics of neural networks ANN, one of the foundational methods in ML, is a network or circuit of artificial nodes that have collective computational properties. Starting from 1950s, McCulloch and Pitts pioneered the subject of ANN by creating a computational model for it[85]. Later on, Rosenblatt invented the concept of perceptron[86], which was described as a statistical model for understanding cognitive systems and became the building block of ANN. Until after one generation, ANN re- gained the attention of scientists because of Hopfield’s work where he introduced non-linearity with coupling between input and output units to an old ANN architecture and therefore endorsed new computational flexibility to this model[87]. Since then, the research interest in ANN has grown steadily and rapidly. Inspired by the biological neurons in a brain, where the electrochemical signal propagates from dendrites (inputs) to axon terminals (outputs) through the myelinated axon (hidden units), a typical neural network has three basic layer-structured components called input layer, hidden layer and output layer, as shown in fig. 2.2. In each layer, there are computing units called neuron. The inter-layer connections represent the flow of data from input to output units. The edges connecting layers are weights and biases to be applied to data emitted from each neuron. Starting from the input layer, a linear transformation is applied to the data, then the neurons in the hidden layer apply a non-linear transformation, often called activation function, to the received data and pass it to the next layer. Note, the number hidden layers and the number of neurons in each layer is not a fixed number, it often depends on the complexity of the problem. As an example, in figure 2.2 we use two, four and three neurons for three layers, respectively. In chapter 4, we will use single neuron for the input, many neurons in the hidden and output layers for our task. 14 Figure 1.2: A schematic of a neural network. The input layer (black) consists of two units, taking in x1 and x2, respectively. The hidden layer (blue) has 4 neurons and the output layer (red) gives an output vector (y1, y2, y3). The blue and red edges connecting layers represent the weights, each node also applies a bias to its input. With the above network configuration and given the input data point _ ⃗ = (_ a ,_ 4 ), the output c⃗ = (c a ,c 4 ,c d ) can be conveniently written as: c⃗ = e a fKe g _ ⃗+W h ⃗ g L+W h ⃗ a (1.6) where e g is a 4-by-2 matrix and W h ⃗ g is a bias vector with length of 4, e a is a 3-by-4 matrix and W h ⃗ a is the corresponding bias for output layer, f:j k →j k is the activation function mentioned above. With this mathematical formulation of the output, one can construct an objective function, i.e. the loss function, which is often defined as certain distance measure between the ANN output and the target output, i.e. true value corresponding to the input data point for a regression task, or the label of the input for a classification task * , then use a proper optimizer to find the weights and biases that minimize the loss function. ANN as an influential tool in multidisciplinary research, the use of it in MD simulations is quite common these days. ANN trained with data from quantum MD calculations are replacing the slow process of force field development[88], some newly developed ANN force fields have been used in MD simulation to study the material synthesis and characterization[32]. ANN is also making an impact on MD data analysis. For example, generative models based on variational * For the classification task, one often needs to do non-linear transformation for the output layer to have predicted label as the final output. Sigmoid or Softmax function are popular choices. Hidden Output Input x 1 x 2 y 1 y 2 y 3 15 autoencoder (VAE), adversarial networks, and normalizing flows have been used to solve inverse problems of materials design from material properties. Two leading exemplars are the VAE models for the design of drug molecules with targeted properties[89] and the study of structural transformation pathways to design novel heterostructures coupling semiconducting (2H) and metallic (1T) phases in strained two-dimensional MoWSe2[90]. Another notable example of ANN in MD is the application of Restricted Boltzmann Machine (RBM) and its enhancement called Limited Boltzmann Machine (LBM) on a quantum annealer to model the synthesis of a MoS2 monolayer by chemical vapor deposition[91]. In this thesis, we report an ANN-based approach to re-formulate MD problem into a boundary value task described by PDEs, and we solve for atomic trajectories in a time-parallel fashion. 1.6 Design of kirigami metamaterials with deep neural network As previously mentioned, TMDs have been highlighted in recent studies for their promising applications in various fields such as catalysis[92], optoelectronic and spintronic devices[93]. Their outstanding properties and the controllability of these characteristics are becoming the spotlight of research in Nanoscience. Researchers have created sensors, high- performance electronic devices using TMDs. For example, Tang et al have grown polypyrrole nanofilms on MoS2 monolayers and fabricated high-performance supercapacitor electrodes. They demonstrated efficient charge storage and transport of both electrons and electrolytes[12]. It has also been shown that the electronic and mechanical properties of these TMD’s can be easily regulated with artificial defects, mechanical straining or by creating hybrid material consisting of vertically stacked heterogenous TMDs. Apte et al. have shown experimentally that under strain MoWSe2 heterostructure goes structural transformation near crack tip from semi-conducting 2H phase to 1T metallic phase, which increases its fracture toughness[18]. 16 Mechanical properties of TMDs can be engineered by creating structural hierarchy inside the material instead of changing its chemical composition[94]. In fact, structural hierarchy has been used to design mechanical metamaterials with negative Poisson’s ratio[95], negligible shear modulus and negative compressibility[96], [97], where truss structures such as kagomé, octet truss is used to synthesize structural hierarchy at atomic level inside the material. Although in general TMDs are brittle in nature, they show high ductility when they are created using the ancient art of krigami design[84]. However, the mechanical properties of these kirigami-based 2D materials are sensitive to the location, geometry and total number of cuts and their interaction among each other[98]. This happens because carefully designed kirigami patterns causes both in-plane and out of plane deformation under strain, which results into high searchability of the materials. It has been shown experimentally that these materials can withstand a strain up to 50%[99]. However, known theoretical results of crack propagation is unable to explain the failure behavior of these 2D kirigami materials since these theoretical results are applicable to single crack inside a material and are unable to explain interactions amongst cuts created in kirigami structures. Further, the search space of the kirigami structures increases exponentially with the configuration of cuts inside the material, which makes the problem of identifying kirigami structures that has maximum searchability an expensive task both experimentally and computationally. In recent years, ML methods have shown tremendous success in material science domain for tasks such as designing data-driven force field models for materials simulation, building ML models to predict materials properly from their atomic structure and even generative models to create atomic structure for materials with a given set of properties. For example, neural-network and gaussian processes has been used to develop data driven interatomic potential for classical MD simulation that can be used to simulate system consisting on millions of atoms at first principal level accuracy[31]. Similarly, ML models has been designed to predict wide range of material 17 properties like band gap, elastic constant, dielectric constants and thermo-electric properties. The input to these models is some representation of the local atomic structure in the material known as atomic fingerprint, and the output is material property. On the other hand, generative models such as VAE and generative adversarial networks (GAN) have shown success in much harder task of creating martials atomic structure for a desired property under consideration, which is also known as inverse design problem[100]. In many situations we are only interested in materials with desired properties such as TMDs with maximum thermoelectric properties, and optimal experimental condition for synthesis or in our current problem: 2D kirigami design with maximum searchability. Further, these problems also involve a sequential decision-making process such as in chemical synthesis the order of reactants, catalyst, temperature and pressure profile are crucial for the final product. Reinforcement learning (RL) is another machine learning domain, which is more suitable for these problems involving an optimization task and a sequential decision making under uncertainty. For examples, RL models haven been used to predict reaction pathways, optimal conditions for chemical reactions[101] and even design molecules with desired chemical properties[102]. In this work we have used RL to design 2D kirigami materials of MoS2 with maximum searchability. We have used classical MD simulation method to create the training data for the RL model, where each MoS2 system consists of up to 6 cuts. Here, the location, length and separation between kirigami cuts are variables in each simulation. Specifically, for the placement of each cut in the material we have 13 different choices based on its length and location. In the RL model, placement of cuts inside the MoS2 structures is posed as a sequential decision-making process, where the goal of RL agent is place up to 6 cuts inside MoS2 structures so as to maximize its total reward at the end which is proportional to the stretchability of the material. Details about our simulation search space and RL model can be found in chapter 5. 18 Chapter 2: Dewetting between MoS2 monolayers 2.1 MD Simulation and validation In the MD simulations reported here, REBO potential is used to describe the interaction between Mo and S atoms including intra-layer bonding and inter-layer L-J interactions[49]. TIP4P/2005 force field is used for modeling H2O molecules. The MoS2 REBO potential accounts for changes in local atomic configurations of atoms, the TIP4P/2005 force field correctly models structural and dynamical properties of bulk and nanoconfined water[103]. The interaction between IPA molecules is described by OPLS-AA force field, which is commonly used for a variety of organic molecules[104]. The interaction between MoS2 membranes and water is modeled by L-J potentials between Mo-O and S-O pairs[50], and interactions between Mo-O, Mo-C, S-O, and S- C pairs in IPA and MoS2 are also described by L-J with different parameters. The force fields for MoS2, H2O and IPA also include long-range Coulomb potential between all charged particles. By modeling MoS2-solvent interface[105], we optimize the force field with experimental data on contact angles for H2O and IPA/H2O droplets on an MoS2 substrate. The force field parameters are listed in table 2.1. Coulomb Interaction Atoms Mo S O(H2O) H(H2O) Charge (e) 1.2 -0.6 -1.1128 0.5564 L-J Interaction Atom Pair Mo-O S-O Mo-H S-H m ij (eV) 4.996E-3 6.929E-3 0.000 0.000 nij (Å) 2.869 3.344 0.000 0.000 Table 2.1: Parameters for Coulomb and L-J interactions between H2O and MoS2 We calculated the interface energy between the solvent and MoS2 from contact-angle simulation and compared the results with experiment. The macroscopic contact angle is estimated using Young’s equation: cos(!)= cos(! r )− s t u v , where w is the free energy correction, ; x is the 19 surface tension of the solvent, and j is the radius of the contact area[50]. To set up the contact angle simulation, MoS2 and the solvent are equilibrated separately for 2 nanoseconds, after which the nanodroplets of different sizes are placed on the MoS2 surface at a height of 3 Å. Fig. 2.1 shows the set up for the contact-angle simulation of H2O nanodroplets. After another 2 nanoseconds of relaxation, the system reaches thermal equilibrium and the density distribution is analyzed in cylindrical coordinates. Fig. 2.2 shows the density profile of a H2O droplet. After taking into account the size effect, we find that the contact angle of H2O droplet on MoS2 surface is 98.9°, close to the experimental value 97.8°[106]. Figure 2.1: MD simulation schematic set up for contact angle simulation. O and H atoms are red and white, respectively; S atoms are yellow, and Mo atoms are not shown in here. (a) (b) Figure 2.2: Density profile of a H2O droplet and the contact angle as a function of droplet size. (a) schematic diagram of the contact-angle measurement. The colored map represents the density distribution in cylindrical coordinates. Z is the height from MoS2 substrate and r is the distance from the central axis of the droplet. ! is the contact angle at the interface of H2O and MoS2. (b) contact angle as a function of radius of the contact area 20 MD simulation setup for dewetting is shown in fig. 2.3. Initially, the system consists of a monolayer of IPA/H2O mixture between two atomically-thin MoS2 membranes of dimensions 100 ´ 100 nm 2 . Water and IPA molecules are distributed randomly at a height of about 3Å above an MoS2 [001] surface and the second MoS2 sheet is placed on top of the solvent at a distance of 3Å from liquid molecules. The IPA concentration is 50% by weight. Periodic boundary conditions are applied parallel to the membranes, i.e., along x and y directions, and equations of motion for atoms are integrated with the velocity-Verlet algorithm using a time step of 1 fs. (a) (b) Figure 2.3: MD simulation set ups for (a) H2O monolayer, and (b) IPA/H2O mixture between MoS2 bilayers. For visualization purposes, MoS2 layers are lifted to expose the liquid layer. Oxygen and hydrogen atoms are red and white, Mo and S atoms are pink and yellow, respectively. IPA molecules are brown. Periodic boundary conditions are applied in the x and y directions of MoS2 bilayers. The size of the system is ~100 ´ 100nm 2 in the x-y plane. We performed DFT calculation for the coupling strength between two MoS2 layers and compared with the MD result based on the REBO forcefield. The DFT cohesive energy is 18.6 meV/Å 2 for an interlayer separation of 3.0 Å, which is the equilibrium separation in our MD simulation. The DFT calculation is based on Perdew’s recent vdW functional SCAN+rVV10, which is good at reproducing interlayer non-bonded interactions[107]. The MD result for the cohesive energy is 15.2 meV/ Å 2 . 2.2 Results and discussion We first examine dewetting in the reference system consisting of a water monolayer between a pair of MoS2 membranes. Fig. 2.4(a) shows the initial configuration; (b), (c) and (d) are three snapshots show the break-up of the monolayer into nanodroplets (red) and small dry 21 patches (white). (MoS2 membranes are not displayed here for the sake of clarity). At the onset of dewetting, we observe small dry patches at random locations in the film. They appear due to thermal fluctuations on the picosecond time scale. After 50 ps, the film breaks up into a network of H2O nanodroplets connected by thin H2O channels. The latter disappear after 200 ps, leaving behind isolated nanodroplets and a much larger fraction of dry patches. The snapshot in Fig. 1(d) shows the result of Rayleigh instability causing the break-up of the monolayer into nanodroplets of diameters ranging between 5 and 20 nm and a few isolated short chains of H2O molecules. The whole process of rupturing of the liquid film takes only 500 ps. Figure 2.4: Dewetting of an H2O film between MoS2 membranes. Red regions represent H2O. (a) Top view of the H2O monolayer at t = 0 ps. Dry patches appear spontaneously and start expanding rapidly. (b) H2O molecules form a network structure at t = 50 ps. (c) At t = 200 ps, the network breaks up into H2O droplets. Some of them coalesce to form larger droplets. (d) H2O dewetting at t = 500 ps. (e) Radius of a dry hole in the film expands linearly with time and the slope gives a dewetting speed of 373 m/s. It is theoretically known that the contact line, i.e., the edge of a dry patch, recedes at a constant speed under specific conditions[108]. To examine whether it is true for dewetting of water between MoS2 membranes at nanoscale, we performed another simulation in which the initial configuration had a dry circular hole of radius 5 nm at the center of the monolayer (see the panel insert in (e)). The initial dry patch was created by removing water molecules in the middle of the monolayer. We monitored the growth of the dry patch as a function of time by tracking the distance between several pairs of points on opposite sides of the dry patch. Fig. 2.4(e) shows that the average separation between those pairs of points on the circular patch increases linearly with time. 22 The speed of the dry patch growth, i.e. the dewetting velocity, estimated from the slope is ~ 373 m/s. This is close to the speed of sound wave in air at room temperature. However, the dewetting velocity estimated from Culick’s law (3 = ( 4t yz {| ) a/4 , where the numerator is the interface energy and the denominator, which is the product of density ~ and thickness Q, gives the surface density of the liquid) is 467 m/s. The discrepancy from our MD result arises from the fact that Culick’s law is based solely on the conversion of interface energy into kinetic energy of liquid and does not take into account the deformation of MoS2 membranes. Our simulation also reveals that rapid dewetting of water is accompanied by significant increases in the temperature and pressure of the entire system. Fig. 2.5(a) shows that the “temperature” (measure of the kinetic energy) of water increases from 300 to 360 K within 40ps while the same temperature increase in MoS2 takes place in 250 ps. Water nanodroplets and MoS2 membranes come to an equilibrium after 300 ps, the temperature of the system remains at 360 K. Fig. 2.5(b) shows that pressure in water due to the deformation increases to 2,700 bars in the first 100 ps of dewetting but does not change subsequently. This increase in pressure does not change the density of water nanodroplet, but produces ripples in MoS2 membranes, and dimples are formed on the membranes after dewetting (see Fig. 2.9). Interestingly, a small cluster of water molecules forms a triangular structure in registry with the MoS2 lattice (Fig. 2.10). Figure 2.5: Temperature and pressure of the system during dewetting. (a) “Temperature” of water rises quickly in the first 50 ps of dewetting. Subsequently, the energy transferred from water to MoS2 increases the temperature of MoS2 membranes. The equilibrium is reached after 300 ps. (b) Shows that pressure increases rapidly as the water monolayer breaks up into nanodroplets. 23 Dewetting behavior of IPA/H2O mixtures is significantly different from that of the pure H2O. Fig. 2.6(a) - (d) are MD snapshots of a mixture consisting of 50% IPA by weight. (a) shows the initial configuration of the system. Thermal fluctuations begin to create small tears in the film. (b) indicates that ruptures have grown into relatively large dry patches (white) after 0.5 ns, and these patches keep on expanding with time as shown in fig. 2.6(c) and (d). Liquid nanodroplets appear between the membranes after 10 ns. They are connected by narrow channels of the mixture, see the snapshot in (d). We also notice that the attraction between H2O and IPA causes water molecules to aggregate around bigger IPA molecules instead of being randomly distributed over the entire wet region. Fig. 2.6(e) shows that the average radius of the dry circular patch increases almost linearly with time, albeit much more slowly than the expansion of a dry patch in an H2O film. The slope of the straight line in fig. 2.6(e) gives an estimate of the dewetting velocity to be around 91 m/s, which is significantly slower than the dewetting velocity in the H2O film. This is due to the fact that IPA molecules are much bigger and hence diffuse much more slowly than H2O molecules. Figure 2.6: Dewetting of an IPA/H2O mixture between MoS2 membranes. (a) Snapshot of the IPA/H2O liquid film t = 0 ps. (b) Snapshot of the mixture taken at t = 0.5 ns shows a percolation network in which IPA and H2O molecules have phase separated. IPA molecules are mostly outside and H2O are inside the network. (c) Snapshot of the IPA/H2O mixture t = 2.0 ns shows an increase in the fraction of dry patches. (d) shows the formation of interconnected nanodroplets after 10 ns. (e) Shows a linear increase in the radius of a circular dry patch (inset). The speed of the dry patch in the mixture is significantly less than that in pure water. 24 We have also monitored the temperature and pressure of IPA/H2O mixtures during dewetting. Fig. 2.7(a) shows how the temperature of an IPA/H2O mixture (1:1 ratio by weight) changes with time. The temperature of the solvent is slightly higher than the temperature of the membranes, indicating that they have not reached equilibrium even after 2 ns. Also note that the temperature increases in the mixture and MoS2 membranes are much less than in the case of pure water. Fig. 2.7(b) shows that the pressure in the mixture due to the deformation of MoS2 membranes also increases much more slowly than in the case of water. The pressure goes up to 1900 bars, which is 30% smaller than the pressure in the previously mentioned water case, indicating a weaker intra-layer coupling. Figure 2.7: Temperature and pressure of the IPA/H2O film during dewetting. The temperature (a) and pressure (b) of IPO/H2O change rapidly in the first 500 ps. Subsequently, the temperature of the mixture and MoS2 membranes change slowly and thermal equilibrium is not established even after 3 ns. In contrast, the equilibrium is reached within 200 ps in the absence of IPA. Another apparent difference between the dewetting of H2O and IPA/H2O monolayers is in the growth rates of dry patches. Fig. 2.8 shows the fraction of total dry areas for pure H2O and an IPA/H2O mixture as a function of time. The fraction of the dry patch in IPA/H2O is lower because wet patches are connected by ligaments in the network structure. The growth rate in H2O dewetting plateaus around 200 ps, whereas dry patches continue to increase slowly during dewetting in the IPA/H2O mixture. 25 Figure 2.8: Time evolution of dry patches in H2O and 50% IPA/H2O mixture during dewetting between MoS2 membranes. Dewetting is much slower in the mixture case because of the slow diffusivity of IPA molecules. The fraction of dry patch in the mixture is lower than that of pure H2O, indicating that the mixture covers more surface area of MoS2 The coupling strength of MoS2 bilayer in the presence of H2O is 7.1 meV/Å 2 , 36.6% smaller than MoS2 membranes without H2O molecules. The cohesive energy in the presence of IPA/H2O mixture is 4.6 meV/Å 2 , which is 34.7% smaller than that of MoS2 with pure H2O. (b) (a) (c) Figure 2.9: Deformation of MoS2 membranes caused by dewetting of the H2O film. (a) Dimples are formed in MoS2 membranes after dewetting. (b) Close-up view of deformation in MoS2. (c) Shows a strip cut from the deformed MoS2 cap. Knowing the bending energy and bending force on the strip, we calculate the bending modulus to be 100±20 N/m. 26 (a) (b) (c) (d) Figure 2.10: Water molecules trapped between MoS2 membranes form two types of structures. (MoS2 membranes are lifted to visualize the structures of water.) (a) Water nanodroplet confined inside an MoS2 cage. (b) Oxygen- Oxygen radial distribution function of H2O in the confined nanodroplet. (c) Confined water molecules form a monolayer of “frozen” H2O in registry with the MoS2 lattice, and (d) shows the triangular lattice structure formed by water molecules. In conclusion, MD simulations reveal distinct dewetting processes in H2O and IPA/H2O monolayers confined between MoS2 bilayers. An H2O monolayer spontaneously ruptures into wet and dry patches with a high velocity (~373 m/s) and wet patches agglomerate to form nanodroplets, which cause local deformations in MoS2 membranes. The entire dewetting process takes about 500 ps, and it is accompanied by significant increases in temperature and pressure of both H2O and MoS2. In contrast, the breakup of an IPA/H2O monolayer is much slower, the temperature and pressure increases are much less than those during the dewetting of a water monolayer. The receding velocity of dry patch is ~91 m/s, and it takes >10 ns before an IPA/H2O film finally transforms into a percolating network of islands connected by narrow channels in which H2O molecules diffuse inside and IPA molecules at the peripheries of the network. The surface area covered by the percolating network of IPA/H2O is much larger than the area covered by H2O nanodroplets. This may explain why mixtures of H2O and IPA are more effective than pure water in liquid phase exfoliation of MoS2. 0 1 2 3 4 5 6 r (Å) 0 1 2 3 g O-O (r) 27 Chapter 3: Origami structure 2D materials 3.1 Experiment details and MD simulation methods Starting from Si wafer spin-coated with polyvinyl alcohol (PVA) layer, we transfer the 2D material onto the wafer, then electron beam irradiation is used to crosslink a selected area of PVA layer to make it insoluble in water. The subsequent immersion into water will remove the unexposed PVA region under the material and make part of the 2D material free-standing. By applying a water flow, the free-moving material will be folded along the edges to a well-defined position and direction, depending on how the electron beam is applied. A schematic illustration is shown in fig. 3.1. Figure 3.1: Demonstration of the folding technique. a) A step-by-step illustration of the folding technique; b) A few-layer graphene flake on PVA/SiO2/Si substrate before and after folding. The AFM image inset shows the folding crease is well aligned with the user-defined PVA edge. c) Deterministic folding of patterned CVD graphene monolayer array under different scales. Dashed and solid arrows indicate the direction of the triangle tip before and after folding, respectively. There is a rectangular PVA layer under each triangular graphene flake, which can be observed by trained eyes. d) An AFM image of a folded graphene sample (after PVA layer is removed) showing the high-quality interface, as well as the well-defined folding edge and folding angle. The inset is the optical microscope image of the same sample, where the area for taking AFM imaging is outlined by a dash line. To mimic the folding process of the material, to understand the mechanism and energetics of the folding, we perform MD simulations of nano-fluid assisted folding of graphene. The initial setup of the MD simulation is shown in fig. 3.2. The system consists of a triangular patch of graphene and 0.27 million water molecules in an MD box of dimensions 24 ´ 24 ´ 12 nm 3 . The 28 length of the equilateral triangular graphene sheet is 20 nm. Carbon atoms in graphene interact via the Adaptive Intermolecular Reactive Empirical Bond Order (AIREBO) potential[109] and we use the transferable intermolecular potential TIP4P/2005 for H2O molecules[103]. The interaction between graphene and H2O is modeled by L-J potential, whose parameters are optimized with experimental data on contact angle for H2O droplets on a graphene substrate[110]. Figure 3.2: MD simulation setup. (a) The simulation box is full of water molecules and in the middle there is a triangular graphene sheet. White spheres are hydrogen, red is oxygen and cyan is carbon atoms. Only half of water is shown in the figure for visualization purpose. We apply constant pressure P = 0.08 GPa on water from right side of the box as the yellow arrow indicates. (b) Top view of graphene sheet in the system. Dashed lines are the boundaries between harmonically-fixed and free-to-move C atoms. The red line is parallel to the left edge of the graphene sheet, the black line is a “tilted” boundary. # is the orientation angle of the boundary with respect to the left edge. In the simulation we have chosen # = 0˚ and 12˚. The hypothesis is that cross-linking of PVA due to irradiation makes the irradiated half of PVA hydrophobic, thereby enhancing the adhesion energy between graphene and irradiated PVA. This enhancement is taken into account by applying a harmonic force to one half of the graphene sample, as shown in fig. 3.2. The spring constant of the harmonic force is taken to be 0.16 nN/Å. If we take spring extension to be 2Å, then the corresponding energy is 0.4 J or 16 kBT. This agrees with experimental estimates of the adhesion energy between graphene and cross-linked PVA. In our simulation we also consider the orientation of the boundary between harmonically fixed and free graphene atoms, the orientation angles are chosen to be # = 0˚ and 12˚, as in fig. 3.2(b). First, an MD simulation is carried out to thermalize the system of graphene immersed in water at 300 K. Periodic boundary conditions are applied in every instance and equations of motion 29 for atoms are integrated with the velocity-Verlet algorithm using a time step of 1 fs. Subsequently, the system is coupled to a heat bath at temperature T = 300 K and the unconstrained part of graphene is subjected to water flow (average speed ~ 100 m/s) at constant pressure P = 0.08 GPa. We have also carried out another simulation with average flow speed ~ 50 m/s and find the same folding phenomena. We do not expect that the water flow speed affects whether folding occurs; it will only affect the rate of the process. We also calculate the free-energy difference, DF, between folded and unfolded states using Jarzynski’s equality[111]: exp(−βW) ÑÑÑÑÑÑÑÑÑÑÑÑÑÑ = exp(−β∆F) (3.1) where W represents the work done during folding of graphene. The free-energy difference between a folded and unfolded (initial) states of graphene is calculated from MD trajectories in the NVT ensemble. 3.2 Results and discussion We first monitored under the fluidic force the increase of the graphene folding angle, which is defined as the angle between current position of the free-standing graphene portion with respect to its initial position. A constant pressure of 0.08 GPa caused the water to flow with a steadily increasing speed. When the front of the flow hits the graphene, water molecules push and carry the graphene patch, driving it to fold, therefore the folding speed is highly correlated with the speed of the water flow. However, a noticeable jump of the folding angle happens at ~150 ps where the folding angle reaches around 160º. It is in fact because of the hydrophobicity of the graphene[112]. When the folding angle is close to 180 º, the sandwiched water molecules get expelled quickly, then the freely-moving graphene portion combines with the fixed graphene portion due to the water pressure difference. This phenomenon can be further confirmed in fig. 3.4(c). 30 Figure 3.3: folding angle of graphene as a function of simulation time. (a) The graphene patch is folded by fluidic force gradually. (b), (c) and (d) are side views of the simulated water (red and yellow) and graphene (cyan) system with a graphene bending angle of 90º, 120º and 180º, respectively. Note, the region colored with yellow is the fraction of water molecules to which we apply a constant pressure to drive the fluidic flow. Fig. 3.4 shows the configurations and energetics of graphene during folding. Snapshots (a)–(d) visualize graphene folded at various angles under the flow of water from the unconstrained side. After folding was complete, we turned off the water flow and found that the folding angles remain unchanged. Fig. 3.4(e) shows the accumulated work on graphene done by water flow, which reaches equilibrium values of 15 and 23 meV for the two graphene configurations, indicating that graphene sheet with boundary orientation angle # = 12° is slightly harder to fold in water than that with # = 0°. Fig. 3.4(f) shows the change of the internal energy of graphene with the increase of the folding angle. The energy barrier of folding a single carbon atom is estimated to be around 2.0 to 3.0 meV, indicating a small out-of-plane bending modulus of graphene. 31 Figure 3.4: Understanding the folding mechanism by molecular dynamics simulation. Snapshots (a), (b), (c), and (d) show graphene bent at angles of 0º, 90º, 150º, and 180º, respectively. Here, oxygen and hydrogen atoms in H2O and carbon atoms in graphene are red, white and cyan spheres, respectively. Yellow spheres represent molecules in the water wake. Note that only half of the water in the MD box is shown in order to visualize the folding of graphene. (e) shows accumulated work per atom done on graphene by water flow as a function of time. Note that φ represents the orientation angle of the boundary between fixed and free portions of graphene with respect to the y axis of the simulation box. (f) shows changes in the internal energy of graphene as a function of the folding angle. The internal energy increases, reaching a maximum value for q = 120º. The potential energy barrier for graphene folding is around 2 to 3 meV. Not only does the water fold the graphene, the graphene also affects the water flow. At the onset of graphene folding, we observe an interesting phenomenon in the formation of a wake in the water jet impinging directly on graphene. The wake evolves into a U-shaped pattern as the folding angle increases. The velocity profile in the wake indicates a Levi distribution and the average wake velocity is 300 m/s. To calculate free energy change of graphene after folding, we performed an ensemble average of Eq. (1) using 20 MD simulations with different initial configurations for each folded state. The free-energy difference per carbon atom increases slowly until the folding angle reaches 90º. Thereafter, the free energy increases more rapidly and reaches a maximum value in the fully folded state. Entropy decrease in graphene sheet after fully folded is 0.6 à â and 0.9 à â for # = 0˚ and 12˚ respectively, indicating that 45% to 60% degrees of freedom for free C atoms has been lost after folding. 32 In conclusion, we developed a fluidic-flow assisted origami technique that can deterministically fold 2D materials and their heterostructures. Using this technique, single-atomic- layer Van der Waals material or their heterostructures with arbitrary scales, thickness, and geometry can be folded without the need of any external supporting polymer layers in the final folded structure. In addition, arrays of patterns can be folded across a large wafer area using this technique and electronic devices that can reconfigure device functions through folding are also demonstrated. The underlying mechanism and energetics of folding are analyzed by molecular dynamics simulations. Such scalable formation of folded Van der Waals material structures with high precision can lead to the creation of new atomic-scale materials and superlattices as well as opening the door to new ways of realizing foldable and reconfigurable electronic devices at micrometer and sub-micrometer scale for memory, sensing, and communication applications. Chapter 4: Nanoindentation on Kirigami MoS2 4.1 MD simulation methods In the nanoindentation simulations reported here, we use a conical diamond indenter (fig. 4.1(c)). To account for the size effect of the indenter, three indenter tips with diameters of 6, 9 and 12 nm are used. The interaction between the indenter and MoS2 is described by LJ potentials between C, Mo and S atoms[113]. However, to avoid the attraction between the indenter and MoS2, in the L-J expression we only keep the 12-th order repulsive components. The interaction parameters including potential depth m HI , characteristic length n and the cutoff distance - Jãå are listed in table 4.1. 33 Figure 4.1: A schematic of nanoindentation simulation (a) Top view of the centrosymmetric nested hexagonal kirigami pattern in a MoS2 monolayer. (b) Top view of square pattern in MoS2. In both (a) and (b), the white region are the places where we have removed molybdenum and sulfur atoms. Yellow region represents the sulfur atoms. (c) Snapshot of the conical diamond indenter we used in the nanoindentation simulation. Cyan spheres are carbon atoms. The indenter angle at the tip is 120 degrees. (d) and (e) are snapshots of conical indenter on top of the MoS2 membrane. Atom Pair Mo-C S-C m ij (meV) 3.325 7.355 nij (Å) 2.818 3.219 -cut (Å) 9.86 11.26 Table 4.1: Parameters for L-J interactions between MoS2 and the indenters We created hexagonal (H) and square (S) kirigami alternating patterns in two MoS2 membranes with a size of 100 ´ 100 nm 2 , as shown in fig. 4.1(a)-(b). The patterns are created by removing specific atoms in the nanosheets based on the predefined patterns. The size of the MoS2 nanosheet is chosen so that it is much larger than the indenter size. After the system is equilibrated at 300K for 2 million MD steps with a timestep of 1 fs, the load is then applied gradually to the basal plane of MoS2. The indenter is pushed down toward the center of the material at a speed of 10 m/s. After every 1 nm the indenter proceeds, it is held fixed in position and the system is relaxed for 100 ps, then the force is calculated and averaged over another 100 ps subsequently. We repeat the process until the MoS2 membrane yields and fractures develop under indentation. Periodic boundary conditions are applied parallel to the basal plane, and a 10-nm wide peripheral region of MoS2 is kept frozen to mimic a rigid platform. Along the z-direction, we use free boundary condion 34 on both sides of the material. The equations of motion for all atoms are integrated using velocity- Verlet algorithm with LAMMPS. 4.2 Results and discussion Fig. 4.2(a) shows the MD results for the load response as a function of nanoindenter displacement, i.e., the indentation depth, for kirigami MoS2 with H patterns with varying size. The patterns shown in panels of fig. 4.2(b) – (d), covers, respectively, 4.1%, 12.7% and 26.1% of the indented area of the material. The load-displacement curves, green, blue and red corresponding to patterns with the same colors, indicate how the mechanical response changes with the area of kirigami patterns. Here, the maximum load is chosen to be just below the critical value at which the membrane fractures. For the same applied load, the indentation depth increases almost linearly with the area of the kirigami pattern. The load-bearing capacity of kirigami structures increases significantly with the area of kirigami patterns. The maximum indentation depth for the smallest pattern (fig. 4.2(b)) is 15 nm and it increases to 30 and 55 nm for kirigami structures shown in fig. 4.2(c) and (d), which indicates that these structures are much more flexible than the pristine MoS2 monolayer. Furthermore, these structures exhibit small hysteresis in the load−displacement curves and suffer less damage than the pristine MoS2 monolayer. Fig. 4.2(a) shows that the hysteresis in the load−displacement curves remains nearly the same for structures (b) and (c). The hysteresis increases for the largest kirigami structure. This trend in hysteresis indicates that pattern (b) and (c) have the same degree of plastic deformation, whereas the largest kirigami structure has undergone more plastic deformation and perhaps other kinds of defects have also formed in the system as a result of much larger indentation depth. 35 Figure 4.2: (a) Load-displacement curves for nanoindentation on MoS2 with 3 kirigami structures of different area fractions. As a reference, the response of the MoS2 monolayer without the kirgami structure is shown in fig. S2 in the supplementary material. The green curve is the indentation response curve of the kirigami pattern shown in (b). The blue and red curves correspond to patterns shown in (c) and (d). The arrows in (a) pointing upward and downward indicate loading and unloading processes, respectively. The hysteresis in load-displacement curves indicates irreversible damage in MoS2 kirigami structures during nanoindentaion. The MoS2 kirigami pattern in (d) can be indented much more with a smaller load than the pristine MoS2 monolayer. Evidently, increasing the area of kirigami pattern has a significant effect on the out-of- plane hardness of MoS2. We have calculated the hardness of these kirigami structures from ç = é èêë í , where ì îïñ is the maximum applied load, P is the effective indented area. Knowing P and ì îïñ , we find that the hardness decreases from 40.9 MPa to 23.6 MPa and 11.5 MPa as the kirigami area fraction increases from 4.1% to 12.7% and 26.0%. These values are much smaller than the hardness 519.6 MPa for a pristine MoS2 monolayer. We have calculated elastic modulii E < for these three structures from E < = √ ò 4 ô Y a √ í , where ö, the slope of the load-displacement curve at the maximum load upon unloading, also represents the stiffness of the MoS2 kirigami structure. The values of S are 6.4 N/m, 3.7 N/m, and 1.8 N/m respectively for kirigami patterns in fig. 4.2 (b), (c) and (d). The quantity ^ accounts for the shape correction factor of the indenter, which is 1.0 for the conical indenter we use in the simulations. The elastic moduli for kirigami structures shown in panels (b) – (d) are 141.2, 81.6 and 39.7 MPa, respectively. Fig. 4.3 displays MD results for load against indenter depth for MoS2 layer with S patterns. The sturctures shown in panels fig. 4.3 (b), (c) and (d) cover respectively 5.4%, 13.3% and 24.7% 36 of the total area of the affected MoS2 surface. Again, we observe that the out-of-plane flexibility of kirgami structures increases with the area of kirigami patterns. The maximum indentation depth before the failure of MoS2 is 16 nm for the pattern shown in fig. 4.3 (c) while the indenter can go as deep as 38 nm for the pattern in (d) before cracks develop at the corners of the pattern. The hardness of the three S kirigami patterns decreases from 42.6 MPa to 10.5 MPa as the kirigami area increases from 5.4% to 24.7%. These values are again more than an order of magnitude smaller than that of a pristine MoS2 monolayer, but similar to the hardness of H patterns shown in fig. 4.2 (b), (c) and (d). The trend in the Young’s modulus for the three S structures confirm the similarity with H patterns. The hardness values are 7.3N/m, 3.4 N/m and 1.8 N/m respectively, for kirigami patterns (b) – (d) in fig. 4.3. The elastic moduli for the three S kirigami structures are 147.6 MPa, 68.8 MPa, and 36.3 MPa, respectively for (b), (c) and (d), very similar to the Young’s moduli of H kirigami structures in fig. 4.2. Figure 4.3: (a) Load-displacement curves for nanoindentation on square MoS2 kirigami patterns shown (b), (c) and (d). The upward and downward arrows in each curve indicate loading and unloading processes, respectively. The area of kirigmai patterns in (b), (c), and (d) are 5.4%, 13.3%, and 24.7%, respectively, of the total area of the MoS2 monolayer. As in the case of hexagonal kirigami patterns, the hysteresis and hardness decrease significantly with an increase in the kirigami area. We have investigated the effect of the indenter size on the indentation response of the largest kirigami structures, where the area of the krigami pattern is ~ 25% of the total area of the MoS2 monolayer for both H and S cases. As shown in fig. 4.4, three conical indenters of diameter 60 Å, 90 Å, and 120 Å were applied in the simulations. The load-displacement curves for H and S 37 patterns are shown in fig. 4.4 (a) and (b), respectively. In the case of the H pattern, the effect of the indenter size is insignificant for indentation depths of 30 Å. Beyond that indentation depth, the force is larger for bigger indenters both in the loading and unloading cycles. The hardness, elastic modulus and stiffness for the three indenter sizes are listed in table 4.2. The load-displacement curves for the S pattern are nearly the same for the three indenters, which implies that the indenter size has an insignificant effect on elastic properties, plastic deformation and defect generation in S kirigami structures. Figure 4.4: Effect of the indenter size on the two MoS2 kirigami structures. (a) Load vs indentation depth for the hexagonal kirigami structure indented with indenters of radii 30 Å, 45 Å and 60 Å. (b) Load vs indentation depth curves for the square kirigami structure indented with the same set of indenters. The area fractions of kirigami patterns in panels (a) and (b) are ~ 25%. As each indenter is retracted, the force on the indenter due to MoS2 deformation reduces to zero slowly. This indicates that MoS2 membrane is in contact with the indenter during unloading. The hysteresis indicates plastic deformation in kirigami structures even after the indenter is completely retracted. The indenter force at each indentation depth is more or less the same for different indenter sizes, which implies that the indenter size has an insignificant effect on the results. The elastic force arises mostly from the deformation of kirigami pattern outside the area of direct contact with the conical indenter. The kirigami structures are extremely flexible and can withstand large strains and stresses without succumbing to fracture. We have calculated stress distribution in kirigami membranes at 38 various indentation depths. Fig. 4.5 shows the stress distribution in H and S kirigami membranes with the same indenter of diameter 60 Å. The area fractions of both membranes are nearly the same (~25%). Fig. 4.5 (a) and (b) show snapshots of the two membranes under the maximum loads. Beyond these loads, the ligaments of kirigami structures fracture and the membranes collapse. At the maximum applied loads, the indentation depth for the H pattern (55 nm) is larger than that of the S (38 nm) structure. Stress distributions in the H and S kirigami at the maximum applied loads are shown in fig. 4.5 (c) and (d). These are atomic-level stresses calculated from the virial expression for the stress tensor. Tensile and compressive stresses are accumulated in the ligaments and corners of both kirigami patterns. The S pattern is evidently under larger compressive stress than the hexagonal kirigami structure due to the different arrangement of the patterns. The stress is negligible outside the indented area. Figure 4.5: Stress distribution in the deformed MoS2 membrane. (a) and (c) are side views of the deformed MoS2 kirigami structures. MoS2 layers form an elastic pyramid-spring under nanoindentation. (Indenter is removed from the view for better visualization of MoS2). Stress is accumulated symmetrically around the corners and arms of kirigami patterns and goes up to 10 GPa. Panels (b) and (d) show top views of deformed MoS2 under nanoindentation. (Stress is calculated using per-atom volume and per-atom stress tensor averaged over 0.2 million MD steps.) Beyond the elastic regime, compressive stresses tend to affect the corners and ligaments more than the rest of hexagonal and square kirigami patterns. Fig. 4.6(a) shows atomic 39 configurations at a corner of a hexagonal kirigami structure before (left panel) and after (middle panel) indentation. Fig. 4.6(b) presents close-up views of a corner in the square kirigami structure before (left panel) and after (middle panel) indentation. As the applied load is increased, stresses weaken Mo-S bonds in the corner regions. Close to the maximum applied load, large compressive stresses break Mo-S bonds in both hexagonal and square kirigami structures. Most of the broken bonds are at the corners and a few of them in the ligaments of the kirgami patterns. These broken bonds give rise to point defects and nanoscale cracks. These defects do not heal upon unloading and nanocracks also remain at the corners of kirigami patterns. We have estimated the density of defects by monitoring the number of broken bonds in the krigami structures both during loading and unloading phases. These results are presented in fig. 4.6(c). Blue curves show how the defect density increases upon loading in the hexagonal kirigami structure. The flat portion of the blue curve indicates that the defect density remains constant in the unloading phase. The lower red curve shows how the defect density increases with an increase in the applied load, and the upper red curve shows that the defect density remains constant during the entire unloading phase. Thus, the indentation leaves permanent damage in the kirigami structures. Figure 4.6: Evolution of defects induced by nanoindentation in MoS2 kirigami structures. (a) close-up views of a corner in a hexagonal kirigami pattern before indentation (left panel) and after indentation (middle panel). There are broken bonds and a nanocrack at the corner after the indenter is removed and the system is relaxed. (b) shows irreversible damage caused by indentation in a square kirgami structure. The left panel shows the atomic configuration at a corner of the square kirgami pattern before indentation. The middle panel is a close-up view of the same corner after indentation. Point defects shown in the middle panel persist even after the indenter has completely retracted and the kirigami structure has been relaxed for 2 million MD steps. The defect density in the two systems is shown in (c). Blue curves represent defect density in the hexagonal kirigami structure in figure 1(d) during loading (lower curve) and unloading (flat upper curve) phases. Red curves in (c) show the defect density during loading and unloading in the square kirigami structure in figure 4.3(d). 40 H Kirigami area fraction Stiffness (N•m -1 ) Hardness (MPa) Elastic modulus (MPa) 0.0% 24.9 160.1 549.2 4.1% 6.4 40.9 141.2 7.9% 4.3 27.5 94.8 12.7% 3.7 23.6 81.6 18.8% 2.9 18.5 64.0 26.0% 1.8 11.5 39.7 Table 4.2: Mechanical properties for H patterns S Kirigami area fraction Stiffness (N•m -1 ) Hardness (MPa) Elastic modulus (MPa) 5.1% 7.3 42.6 147.2 13.3% 3.4 20.7 68.6 24.7% 1.8 10.5 36.3 Table 4.3: Mechanical properties for S patterns Chapter 5: Kirigami material design with deep neural networks 5.1 MD simulation and machine learning formulation In our simulation, MoS2 krigami structures is created by removing Mo and S atoms from the MoS2 nanosheet based on a predefined pattern and such that stoichiometry of the system is maintained. Each cut has the same width of 5 nm and each nanosheet has dimension 20 ´ 30 nm 2 consisting of ~20,000 atoms. The created systems are first relaxed to a minimum energy configuration using conjugate gradient and then thermalized at 100K for 200 ps. The system is then subjected to uniaxial tensile deformation which consists of a loading phase and a relaxation phase. During the loading phase, system is homogenously expanded along y-axis for 1000 ps, which is followed by relaxation phase where the entire system is relaxed for 5 ps, and after that the whole process is repeated again. Periodic boundary condition is applied only along the y-axis, and the tensile stress of the material is computed after the relaxation phase. The simulations are performed using Stillinger-Weber potential[47], [114] with LAMMPS. 41 We compute the stress of the material using symmetric kinetic energy tensor and the virial tensor, which can be expressed as: ì SY = ∑ + H 3 H,S 3 H,Y ú Hùg û + ∑ - I,S ) I,Y ú Iùg û (5.1) where ì SY and û are the stress tensor components and volume of the simulation box, respectively. The summation runs over all the atoms, 3 H,S and 3 H,Y are the components of the velocity for i th atom, - I,S and ) I,Y are the components of the coordinates and force components applied to j th atom. A schematic illustration of the MoS2 kirigami structure, and its corresponding stress-strain curve is shown in fig. 5.1(a) and (b). Here, the stretchability of the material is defined as the maximum strain at which material fails - the point of sudden drop in stress in the stress-strain curve of the material. We can observe from fig. 5.1(b) that failure point of these materials is very sensitive to the topology of the cut patterns inside the material. In fact, their stretchability can be shown to be the function of total number cuts, cut length and their location on the nanosheet. Fig. 5.1(c) shows some possible choices of kirigami cut patterns in a 20 ´ 30 nm 2 MoS2 nanosheet. The strategy for creating cut is that the system is first vertically partitioned the nanosheet into n rows of equal width, where each of these rows may contain a horizontal cut. In these rows, a single cut can be placed at one of the 4 or 6 equally-spaced points chosen horizontally along the central axis of the row (fig. 5.1(a)), and the length of these inserted cut is either 5, 10 or 15 nm. Further, if the end point of the inserted cut extends outside of the MoS2 nanosheet width (20 nm), it is wrapped back from the other side. Thus, we have 13 different choices to place a cut in each row: 4 different location ´ 3 cut length + no cut, which makes the total number of possible krigami structures 13 n , where n is the total number of vertical rows in MoS2 nanosheet. We can observe that state space of structure grows exponentially with the number of rows containing cuts, which for n = 4 is 28,561 and for n = 6 it becomes 4,826,809. Due to the vast search space, it is not feasible to find structures 42 with high stretchability via experiments or MD simulations. Instead of brute-force search to find optimal structures, we have created a RL model that propose structure with high stretchability after training. The data for the RL model is created by doing uniaxial tensile MD simulations of all the 28,561 structures for n = 4 and randomly sampled 38,406 structures for n = 6. The details of the MD simulation are given below. Figure 5.1: illustration of MoS2 kirigami structures and their corresponding stress-strain curves. (a) the horizontal lines divide the MoS2 nanosheet into independent 6 rows, in each row we define 13 different configuration of a cut slit. All cuts in the 6 rows combined form a kirigami pattern in the MoS2 nanosheet. (b) The stress-strain curve corresponding to the kirigami structure in (a). (c) examples of possible kirigami patterns in our design space. Our RL model consists of an agent and a reward model, and together they propose kirigami structures as a sequential decision-making process. Here, synthesis of a kirigami structure from a pristine MoS2 nanosheet is called an episode of RL, where the length of each episode is kept fixed at T = 4 or 6 (equal to the total number of vertical rows). During an episode at each time step t, the RL agent takes a decision about the location and length of next cut in the t th row from the 13 available choices (action space), where the input to the RL agent is the structure from the (/−1) åü timestep. Using the RL agent proposed action, a cut is created at the specified location and length in the / åü row, which serves as input state to the RL agent at next time step. Further at each time step the model gives a reward to the RL agent, which at intermediate time steps /<° is 0 and at terminal state (/= °) is proportional to the stretchability of the material following: 43 jQe,-¢ (-) = £ 0 ;/≤°−1 ¶- /= ° ,ß¢ ® = 15% ∝0.1ö ; /= ° ,ß¢ ö >15%,eℎQ-Q ® = stretchability (5.2) Fig. 5.2(a) and (d) show a schematic of our RL agent, reward model and structure design in each episode of RL, respectively. The reward model is a deep convolutional neural network (CNN), where the input is a 1´16´16 binary tensor representing the kirigami structures, and the model’s output is stretchability. The RL agent is a Deep Q-Network (DQN), whose input is a 2´16´16 tensor that summarizes the intermediate state (® å ) of the structure at timestep t, and output the state-action value function, N(® å ,, å ) for all the actions we can take from ® å , which is 13 different cut chooses to place in the material in the / åü row. Here, N(® å ,, å ) is an estimate of the total expected reward for each state-action pair from / åü time-step to the end of episode under a given policy, and the goal of the DQN is to learn a policy (∂) that maximizes the N(® å ,, å ) for each state-action pair, see equation 5.2. The details about training of DQN and reward model is given in table 5.1. N(® å ,, å )= max ò E[- å +; - å∏a +; 4 - å∏4 +⋯; ∫RåRa - ú ]≈ max ò E º Ωæ- å +; max ï ø¿¡ N(® å∏a ,, å∏a )¬ (5.3) Figure 5.2: An illustration of the RL model. (a) a schematic of the RL agent, the input to the agent is a binary tensor representing the kirigami structures and the output is the corresponding stretchability, i.e. the maximum strain that this kirigami MoS2 can withstand. (b) loss function against number of training epochs. (c) the expected reward as a function of training epochs. 44 Initialize Replay memory of size N with by random samples of tuples (® å ,, å ,- å ,_ å∏a ). Create target DQN, N(! √ ) , network with random weight initialization, ! √ Create policy DQN, N(! ∗ ), network with random weight initialization, ! ∗ ≈¶- Q∆[®¶¢Q = 1,« Initialize sequence ® g = _ a , ℎQ-Q _ a [® «¶ö 4 ®/-»…/»-Q e[/ℎ ß¶ …»/ ≈¶- /[+Q /= 1,° - = -,ß¢¶+ ß»+WQ- WQ/eQQß (0,1) [) - <Q∆®[ [¶ß (m ): , å = …ℎ¶¶®Q , ,…/[¶ß -,ß¢¶+ c )-¶+ P Q ®Q: , å = ,-f+, _ ï∈í N(® å ,, å ;! ∗ ) )-¶+ ∆¶ […c ßQ/e¶-à E_Q…»/Q ,…/[¶ß , å ,¶W®Q-3Q -Qe,-¢ - å ,ß¢ ßQ_/ ®/,/Q ® å∏a (…»/ [ß / åü -¶e) ,¢¢ /»∆ Q (® å ,, å ,- å ,® å∏a ) to the reply memory ö,+∆ Q , +[ß[W,/…ℎ )-¶+ /ℎQ -Q∆ c +Q+¶-c c I = Ã - I [) \ = ° - I +; max ï N(® Õ∏4 ,, ,! √ ) [) \ ≠° Minimize the loss between target and policy network- œc I −NK® I ,, I ;! ∗ L– 4 – by one step of gradient decent. —) +¶¢(Q∆[®¶¢Q,8)== 0: 8∆¢,/Q /,-fQ- ßQ/e¶-à: N(! ∗ )= N(! √ ) Table 5.1: DQN training pseudocode with experienced replay 5.2 Results and discussion Two RL agent are constructed using the same DQN architecture for n = 4 and 6, respectively such that they are trained to learn a policy that proposes structures with at least S = 15% called threshold stretchability ö ∫ . Loss and the expected reward per episode during training for these models are given in fig. 5.2(b) and (c) respectively. After training, RL agents follow an “-greedy policy to sample new structures in each episode, where at each timestep it chooses an action that has maximum Q(s,a) with 1−“ probability or a random action with “ probability. Fig. 5.3(a) shows the true distribution of stretchability of 28,561 structures for n = 4, and a histogram of stretchability of 500 structures sampled using a random search and by the RL agent for n = 4, in which case each structure is chosen with equal probability. We can observe from fig. 45 5.3(a) that the distribution of ö is highly sparse, where 90% of structures have ö ≤15%. Thus, a random search is not able to find structures with superior mechanical properties, which is evident since majority of the structures have ö ≤15%. Further, the random search is not able to find any structure closer to maximum S (ö îïñ ). On the other hand, after training the RL agent proposes structures with high S 90% of the time. In fact, most of the structures proposed by the RL agent have ö ≥20%, which is closer to the true distribution structures nearer ö îïñ for n = 4. Fig. 5.3(b) and (c) also show some structures proposed by the RL agents and randomly selected along with their S values. Figure 5.3: distribution of max strain and comparison between RL-generated and randomly selected structures. (a) shows the distribution of the maximum strains for the entire state space for n = 4 and sampled space for n = 6. For both cases, the majority of the structures have max strain less than 20%. (b) distribution of the max strain of 500 RL generated kirigami structures, along with a visualization of some of the structures. (c) distribution of 500 kirigami structures uniformly sampled from the database, and some typical structures with their strechability. The sequence in which we place cut inside the material affects its S heavily. In fact, we can quantify the importance of each cut in terms of its state-action value function, N(® å ,, å ). Fig. 5.4(a) shows the principal component analysis (PCA) of the learnt feature in the last layer of DQN for all the intermediate structures and the action associated with it that has maximum N(® å ,, å ) at time /= 2 and 3 for ß = 6, that means for the second and third cut places inside the material. We can observe that there are very few design choices available with respect to the placement of the 2 nd and 3 rd cut inside the material, which will eventually lead to the creation of structure with high S 46 values. Two such intermediate structures that high S value is also shown in fig. 5(a) and the creation of kirigami structure by the RL agent via these two intermediate structures is shown in fig. 5(b). Using N(® å ,, å ) we can also determine the best location and length to place the very first cut inside the material. Figure 5.4: (a) Principle component of the state-action value function against cut sequence. (b) illustration of the creation of kirigmai structure by the RL agent. In the case of n = 6, due to the large search space consisting of millions of structures, we have only used 70,000 structures (1.45% of the state space) to train the RL model, and after training used the RL agent to find structures with highest S as shown in fig. 5.3(b). We have further validated the S values of these structures using uniaxial tensile MD simulation, and their true and RL model predicted S values agree. These structures show ductile fracture and thus have high S, whereas structures with smaller ö ≤10% shows brittle failure. Fig. 5.5 (b), (c), (e), (f) and (g) show the atomistic stress distribution for two stretched kirigami MoS2 nanosheets and the stress- strain curves for a system with lower and higher S, respectively. In case of structures with low S, e.g. fig. 5.5(a)-(c), the local tensile stress concentrated near cut edges and connection bridges increases with applied strain. When any of these local stress increases above a critical value (bond strength of Mo-S bond), it results into Mo-S bond breaking and rapid crack propagation originating from those high-stress regions (red region in fig. 5.5(b) and (c)). The cracks also result into sudden drop in the system stress as in (g) and shows brittle fracture. However, in the case of systems with 47 high S value (fig. 5.5 (d)-(f)), under strain beside local tensile stress conc. near cut edges, we also observe the presence of compressive stresses in the intermediate regions of the cut both along x and y direction (fig. 5.5(e) and (f)). Due the absence of constraints along z axis, these compressive stress causes out of plane motion (fig. 5.5(d)) along the z-axis in those regions, which results into releasing system’s stored energy and delays the tensile stress near cut edges to reach Mo-S bond strength, and we observe high S value. Further, these out-of-plane deformation of localized regions inside the systems due to the compressive stress varies depending upon cut design created on the system, which then directly affects its S value (fig. 5.5(g)). Figure 5.5: atomic stress distribution and the stress-strain cuve. (a) sideview of a brittle MoS2 kirigami structure that is about to fracture, no out-of-plane deformation is observed. (b) and (c) different components of the atomic stress, tensile stress accumulates at the weak points of the structure. (d) side view of a ductile MoS2 structure, where out-of-plane deformation happened. (e) and (f) not only tensile stress is concentrating at the corners and bridges of the cut pattern, we also observe compressive stress where out-of-plane bending is observed. (g) the stress-strain curve of the two structures shown in (a) and (d) In conclusion, we observe that RL can be used to generate structures with high stretchability from an extremely large search space consisting of millions of structures, which is otherwise impossible to determine from random search due to the fact very few structures have high stretchability value (< 5%). Further, RL agent also provides insight the important of placement and length of cut inside the material that will eventually create structure with superior mechanical properties in terms of its state-action value function, N(® å ,, å ). 48 Chapter 6: Principle of least action for neural network 6.1 Action-derived MD with neural network Least action is the foundational principle of physics. It pervades classical mechanics[115], special and general theory of relativity[116], quantum mechanics[117] and thermodynamics[118]. Its variant, the principle of least time, unifies ray and wave optics[119]. An artificial neural network (ANN) is a network or circuit of artificial neurons or nodes that have collective computational properties[87]. In this paper, we describe how an ANN learns the principle of least action by minimizing an action: P = ’ D å ¡ å ÷ ¢/ (6.1) and provides the time evolution or dynamics of a many-body system between initial and final configurations at time t0 and t1, respectively. In classical mechanics, the L in the action integral is the Lagrangian, it is the difference between the instantaneous kinetic and potential energies of a system and the classical action can be viewed as the difference between the time-averaged kinetic and potential energies. Action A is the key quantity in the loss function of our neural network. The minimization of the loss function provides the trajectory q(t) and its time derivative ◊̇(/) which collectively represent positions and velocities, respectively, of particles in a system. The principle of least of action embodies a boundary-value problem in which A is minimized with respect to q(t) subject to constraints that the initial and final configurations, q(t0) and q(t1), are fixed[120]. Another formulation arising from the principle of least action is Euler- Lagrange equations, which constitute an initial-value problem. These ordinary differential equations are at the core of MD, the preeminent dynamic simulation approach in physics, materials science, and biology[121]. The essential input to MD is interatomic potential energy from which forces are calculated and Newton’s equations or Hamilton’s equations of motion are integrated 49 over discretized time with a finite-difference scheme[34]. The output is phase-space trajectories {◊(/),◊̇(/)} from which structural, mechanical and thermodynamic properties of the system can be computed and compared with experimental measurements. While the neural-network approach to the principle of least action is new, the use of machine learning (ML) in MD simulations is quite common these days especially in the development of MD force fields[88]. ANN trained by data from quantum mechanical calculations are rapidly replacing the painstakingly slow process of force-field development. ANN-based force fields are now widely used in MD simulation studies of material synthesis and characterization. ML is also making a huge impact on data analytics related to MD simulations. Predictive models based on kernel ridge regression, support vector machines, and random forests have been developed to predict material properties such as band gaps[122], elastic constants[123], thermoelectric properties and so on[124], [125]. ML combined with classical or quantum MD simulations have been used in a variety of other ways in condensed matter physics, chemistry and materials science. For example, Bayesian optimization methods together with MD have been used to discover optimal layered materials for targeted material properties such as electronic band structure and thermal-transport coefficients[43]. Generative models based on the variational autoencoder (VAE), adversarial networks, and normalizing flows have been used to solve inverse problems of materials design from material properties. Two leading exemplars are the VAE models for the design of drug molecules with targeted properties[89] and the study of structural transformation pathways to design novel heterostructures coupling semiconducting (2H) and metallic (1T) phases in strained two-dimensional MoWSe2[90]. Another notable example of ML in MD is the application of Restricted Boltzmann Machine (RBM) and its enhancement called Limited Boltzmann Machine (LBM) on a quantum annealer to model the synthesis of a MoS2 monolayer by chemical vapor deposition[91]. 50 In this paper, we have taken the first step in demonstrating that a feed-forward neural network (NN) can learn the principle of least action and provide atomic trajectories that are in excellent agreement with trajectories obtained by the conventional MD approach based on a time- stepping algorithm for sequential integration of Newton’s equations of motion. We have chosen to train our NN on the Onsager-Machlup (OM) action: ö Ÿ⁄ = ’ ¤G‹+ H ◊ ⃗ ̈ fi (/)+ flû(N(/)) fl◊ ⃗ fi ‹ 4 ú Hùg ‡¢/ ∫ g (6.2) where mi is the mass and ◊ ⃗ ̈ fi (/) is the second-order time derivative (i.e. acceleration) of the position vector ◊ ⃗ H (/) of i th particle. In Eq. (6.2), N(/) collectively represents the coordinates of all the particles in the system at time t and û(N(/)) is the potential energy of the system. This action was derived for irreversible stochastic thermodynamic processes[126]. The genesis of the OM action was Onsager’s seminal papers in which he developed reciprocal relations for transport coefficients and proposed that the probability for fluctuation histories is exponentially small with the exponent proportional to the time integral of dissipation function of the history[127], [128]. Later on, Onsager and Machlup showed that the most probable history is given by the principle of least dissipation in which the dissipation function takes the form of a Lagrangian in classical mechanics and the minimization of eq. (6.2) would produce the most probable trajectory[126]. Here we show how a neural network minimizes ö Ÿ⁄ to produce the entire trajectory of a classical many-body system in one fell swoop. Our NN for the OM action can be easily implemented in any statistical ensemble and it is straightforward to include constraints, invariances, and conservation laws in the neural network. An intriguing aspect of using a neural-network approach to derive the dynamics of a system from an action such as eq. (6.2) is the general possibility of splicing the trajectory into time segments and performing parallel calculations in time. This might be a way to address the long- 51 standing problem in atomistic simulations that they cannot reach long timescales except for a few particular cases where the temporal locality of natural processes and transition-state theory allow the reformulation of the sequential long-time dynamics as a parallel search for low activation- barrier transition events. However, in general, it has been a challenge to simulate events involving slow atomic processes or slow conformational changes such as protein folding[129] or dynamics of recrystallization in phase-change materials[130]. Space-time parallelization through domain decomposition combined with highly efficient algorithms and massively parallel machines would enable large-scale atomistic simulations of complex materials to reach long timescales[42]. 6.2 Neural network structure and simulation methods We train a neural network (NN) to learn the OM action for a many-body system of L-J atoms[131]. The OM action has become a popular optimization approach to solving atomistic boundary-value problems because it is guaranteed to have a global minimum[132]. However, selecting a starting trajectory for the OM action is a non-trivial task. One approach is to start with a straight-line trajectory between the initial and final configurations, express the trajectory in terms of a Fourier series and vary the Fourier coefficients to minimize the action in eq. (6.2) under the constraint of energy conservation. Unsurprisingly, the final trajectory is found to depend heavily upon the initial values of Fourier coefficients[133]. Lee et al. attempted to address this problem by introducing in the cost function an additional term that makes the average kinetic energy consistent with the temperature of the system[134]. Later on, they used conformational space annealing to obtain many low OM action trajectories[135]. Fig. 6.1 shows our feedforward neural network which consists of an input node, a hidden layer with n units, and an output layer with d´N units, where N is the number of atoms. The neural network outputs the trajectories of all the atoms in the system, i.e. cartesian components {◊ H,ñ , ◊ H,Õ , 52 ◊ H,· } for i = 1 to N. We train the neural network to compute configuration-space trajectories of L- J systems consisting of N = 500 atoms in three-dimensional (3D) case. Figure 6.1: A schematic of the neural network. The system configuration for all spatial degrees of freedom is calculated by inserting an instant of time t into the input layer. A single hidden layer (blue) is used and the output layer (red) gives the atomic coordinates whose first and second derivatives provide velocity and acceleration, respectively. The blue and red lines connecting layers represent the weights and each node also has a bias. We evaluate the OM action numerically by discretizing the time integral in eq. (6.2) into ‚ å grid points (eq. (6.4)). The time increment is ∆/= °/‚ å , where T = 1 ps and ‚ å = 25. The loss function L for the neural network includes not only the OM action but also constraints due to energy conservation and boundary conditions: „({e,W }, ‰ a , ‰ 4 , ‰ d , ‰ k ) = ‰ a ö Ÿ⁄ ÂHº + ‰ 4 (N(0)−N g ) 4 + ‰ d (N(°)−N ∫ ) 4 + ‰ k GÊEKN(/),N ̇(/) L−E g Á 4 ∫ åùg (6.3) ö Ÿ⁄ ÂHº = GË(N ̈(/)+ ¢û(N(/)) ¢N )Ë 4 ∆/ ∫ åùg (6.4) EKN(/),N ̇(/) L= G 1 2 + H ◊ ⃗ ̇ H 4 ú Hùa + 1 2 GG4“Í n a4 Î◊ ⃗ H −◊ ⃗ I Î a4 − n Ï Î◊ ⃗ H −◊ ⃗ I Î Ï Ì ú IZH ú Hùa (6.5) where {e,W } represents all the weights and biases in the neural network. For simplicity, we use N(/) to denote N({e,W },/). N g and N ∫ collectively represent initial and final coordinates of all the atoms in the system and ‰’s are the hyper-parameters. Note that the last term in eq. (6.3) is the energy-conservation constraint, where E g is the initial total energy of the system. In eq. (6.5), 53 EKN(/),N ̇(/) L is the total energy with “ and n as parameters of the L-J potential. The parameters l are chosen with care to ensure that all the boundary conditions are well satisfied by the neural network. We find that setting ‰ a = ‰ k = 1 +,_ å (|û(N(/))|) Ô ensures that the neural network sets the initial conditions correctly before minimizing the action. In this work, we use the reduced L-J units, the interaction parameters for the L-J potential are given in table 6.1. Molecular dynamics simulation was initiated from the face-centered cubic (FCC) structure. The reduced number density of the system is 0.787 and the reduced time unit w = 2.068 ps. Periodic boundary conditions were imposed, and equations of motion were integrated with the velocity-Verlet algorithm using a timestep of 2 fs. The system was equilibrated for 2 ns at a reduced temperature of 0.695. We use the coordinates and velocities of atoms in the equilibrated state as the initial configuration as the input to NN, and we run another MD simulation for 0.5 w using the same starting state obtain a final state for our boundary-value problem. Simulation Parameter Value ϵ/à F (Ú) 125.7 σ (nm) 0.3345 m (amu) 39.948 Lattice Constant (nm) 0.575 Table 6.1: L-J interaction parameters and lattice constant In order to obtain an initial guess for the NN output, we need to pre-train our network. The pre-training strategy we use is as follows: we first randomize the parameters {e,W }, then with the knowledge of the initial and final state of the system, we pre-train the network using the following loss function: 54 D ıT| ({e,W })= KN h ⃗ ({e,W },0)−◊ ⃗ g L 4 +KN h ⃗ ({e,W },°)−◊ ⃗ ∫ L 4 +G(N h ⃗ ̇ ({e,W },[)−3 ⃗ g ) 4 ∫ 4 ⁄ Hùg + G (N h ⃗ ̇ ({e,W },\)−3 ⃗ ∫ ) 4 ∫ Iù∫ 4 ⁄ ∏˜å (S6.1) where ◊ ⃗ g and ◊ ⃗ ∫ collectively represent the initial and final positions of atoms, and 3 ⃗ g and 3 ⃗ ∫ are initial and final velocities of atoms in the system, respectively. Note that the first two terms are constraints on initial and final configurations of the system, while the last two terms provide an estimate of trajectories by making particles move at different but constant velocities for the first and second half of the trajectories. During the minimization of the loss, we use Nestrov and Adam optimizer (NADAM)[136], and the gradient of the loss function with respect to parameters is calculated numerically using Google’s JAX library[137]. After training the neural network, we compute atomic trajectories and compare them with the “ground-truth” conventional MD simulations which are performed with the same set of initial conditions. 6.3 Results snd discussion Our NN learns to compute atomic trajectories of 3D L-J systems from the OM least action principle. Here we present results for the system containing 500 atoms, which was trained on a neural network with 125 neurons in the hidden layer for 200,000 epochs. Note, the results in this section are generated on a time domain 0 to T = 0.5t with Nt = 25 time increments which correspond to a time step (~ 40 fs) that is 20 times the MD timestep (2 fs). Fig. 6.2 shows a comparison of four randomly chosen single-particle trajectories in the system over a time period of 1 ps. These trajectories cover a time scale beyond the straight-line ballistic motion and possess sharp changes in direction. It is remarkable that the neural network (NN) is able to reproduce the dynamics at the individual particle level. From this level of matching 55 the MD results for particle trajectories, it is not surprising that the NN results for all other quantities are in excellent agreement with the “ground truth” MD results. Figure 6.2: Four selected NN (solid blue) and MD (solid red) trajectories show that the NN can reproduce non- trivial paths. These are for a time period of 1 ps. In fig 6.3, we compare the results from the NN and MD for some standard simulation quantities. The total energy is part of the loss function, thus we expect it to hold quite well, as shown in fig. 6.3(a), the NN kinetic and potential terms also matches the MD values. The structure of the liquid is characterized by the radial distribution function[138]. g(r) shown in fig. 6.3(b). The NN matches the MD result. Slightly below r = 1, g(r) = 0 because the particles cannot overlap. The NN trajectory maintains this feature as well. The peaks in g(r) corresponding to the minimum energy separation for nearest and second nearest neighbors indicate that the NN is able to obtain the correct liquid structure. To characterize the particle dynamics, we calculated the velocity autocorrelation function[139] (VAF) shown in fig. 6.3(c). The match between NN and MD is good, and the correlation time tc, when the VAF = 0 for the first time, is about 0.16 tc. Beyond tc the velocities and, in general, the particle dynamics are not correlated with the initial dynamics of the particles. This quantity is important in calculating time averages, as a time series of data with points separated by at least tc are not correlated and contribute independently to the average quantity. 56 Figure 6.3: Comparison between the neural network (NN) and MD simulation. (a) Kinetic (blue), potential (green) and total (red) energy per particle as a function of time for MD (square dots) and NN (solid lines). Radial distribution function and velocity autocorrelation function for MD (blue) and NN (red) are shown in (b) and (c), respectively. Fig. 6.4 shows how various contributions to the loss function decrease with the number of epochs. Since the system is in the microcanonical (NVE) ensemble, the constraint of energy conservation is part of the total loss function. All three terms in the loss function decrease rapidly in the first few thousand epochs. The contribution from the OM action is the largest and determines the minimal value of the total loss. The loss function corresponding to boundary-condition constraints decreases rapidly and reaches below 10 -6 . It is not surprising that this loss is so low because these conditions are only on two points in the trajectory. Contrast this with the loss for the energy-conservation constraint, which goes below 10 -5 after 200,000 epochs since the constraint must be satisfied for the whole trajectory. We also show in fig. 6.4 the mean-squared error (MSE) of the position (Q) and momentum (P) with respect to the corresponding MD variables. These errors are computed from trajectories generated by the neural network and the “ground-truth” MD simulations. The general shape of the MSE is similar for the two quantities although they differ by orders of magnitude. We offer two explanations for the rising error per time derivative. First, referring to equations (6.3) - (6.5), the time derivatives are described with fewer parameters than the positions as the biases outside the activation function are killed by the derivative. Second, all quantities in the loss function depend explicitly on position, two depend on the velocity and only one depends on the acceleration. Since we have a finite number of parameters, they can only approximate the function up to some error. 57 These errors magnify as we take higher derivatives. The mean square error for the positions is very small, which as we saw in fig. 6.2, results in excellent matching of the neural network and MD trajectories. The behavior of the momentum as a function of epoch is similar to the position’s, but with a larger MSE at convergence. Figure 6.4: (a) shows how the three components of the loss function decrease with the number of epochs. The boundary condition term (red) and energy conservation term (blue) fall more rapidly than the OM action term (green). (b) Shows the mean square error (MSE) for trajectory (red) and momentum (blue) which was calculated by comparing the MD trajectory to the neural network predicted trajectory. (c) shows differences between MD and NN trajectories (red) and velocities (blue) as a function of time. Here t = 2ps. We also studied how the NN scales with the number of units in the hidden layer and discretization of the time domain on Intel Xeon(R) CPU E5-2640 machine with 2.60 GHz clock speed. Surprisingly, with the help of Google JAX, fig. 6.5(a) shows that increasing the number of neurons in the hidden layer doesn’t affect the run time appreciably. The runtime remains around 210 seconds for 1,000 epochs of training. However, changing the grid size of the time domain increases the computing time significantly, as shown in fig. 6.5(b). Fig. 5(c) shows how the run time scales with the system size. The behavior appears quadratic because the run time is dominated by the force calculation. It’s straightforward to reduce the scaling of the run time with the system size to O(N) by using the linked-list method[131]. 58 Figure 6.5: Scaling of NN for 108-particle system. (a), (b) and (c) are the run time for every 1000 training epoch as a function of the number of units in the hidden layer, the size of the time grid, and the number of atoms in the system, respectively. In this work, we have shown that a feedforward NN can correctly simulate the phase-space movement of a many-body L-J system by minimizing the Onsager-Machlup action. Besides the atomic trajectories, the potential and kinetic energies, radial distribution and velocity autocorrelation functions are in agreement with the corresponding results from the “ground-truth” MD simulation. Our neural network approach is generally applicable to systems with more complex interatomic interactions than the presented L-J potential. Although the time efficiency of our NN model scales as N 2 , we have made the first-step effort to demonstrate the viability of this method. The proposed neural-network approach to principle of least action has a few advantages over the conventional MD method: our NN provides the entire trajectory in one fell swoop over a time period of 1 ps; the time step in the NN approach to action is a factor of 20 larger than the time step in conventional MD simulation; even with a much larger time step, the energy conservation over the entire time domain is as good as in the conventional MD simulation. Going further, the neural-network approach presented here has the potential to be a game changer in more than one way. For example, Hills et al. have developed an algorithm based on the least action principle to predict the dynamics of physical systems using observed data[140]. Our approach combined with theirs can model observational data for different types of physical systems. This can be very useful in predicting the future dynamics of physical systems outside the 59 range of measurements and quite conceivably the two neural networks could help scientists discover new phenomena hidden in data. Another important area in which our neural-network approach can be impactful is imputation of missing data. Suppose we have disconnected trajectories of a physical system and we want to find the missing piece between the end point of the first and the starting point of the second trajectory. Using the algorithm of Hills et al., one can find the Lagrangian of the underlying system. Subsequently, one can use our neural network to interpolate the end point of the first and the starting point of the second trajectory. Chapter 7: Summary and future directions 7.1 Summary 7.1.1 Origami and kirigami based 2D materials In this thesis, we have explored folding and structure engineering of 2D materials via the simple art of origami and kirigami. For the folding of graphene and other 2D materials, we developed a fluidic-flow assisted origami technique that can deterministically fold 2D materials and fabricate heterostructures. We also analyzed the underlying mechanism and energetics of the process via MD simulations. Using this folding method, single-atomic-layer 2D materials with arbitrary geometry and scale can be folded without the need of any supporting layers in the final structure. Such scalable formation of folded 2D structures with high precision can lead to the creation of new atomic-scale materials and superlattices as well as opening the door to new ways of realizing foldable and reconfigurable electronic devices at micrometer and sub-micrometer scale for a wide range of applications. 60 To measure the effect of kirigami on the mechanical properties of 2D materials, we performed MD simulations to study the deformation of kirigami MoS2 monolayer under indentation. We found that the kirigami patterns can dramatically modulate the rigidity of an atomically thin 2D material. A pattern with ~20% area of the MoS2 membrane can reduce the out- of-plane stiffness by 2 order of magnitudes. The indenter size effect is insignificant in the out-of- plane deformation since most of the distortion appears at ligaments and corners of the kirigami patterns that are not in direct contact with the indenter. Via simulation, we identified two types of point defects in the deformed kirigami structure, which prevented the elastic recovery of the MoS2 membrane from the indentation experiment and introduced hysteresis in the load-displacement curve for all of the kirigami MoS2 structures. We also propose nanoindentation experiment and potential applications with this novel kirigami structures in the monolayer MoS2 membrane. 7.1.2 Neural network meets molecular dynamics In this thesis, we have shown that a multi-output feedforward neural network can accurately predict atomic trajectories of a L-J system by optimizing the objective function incorporating Onsager-Machlup action and necessary constraints. In addition to the phase-space atomic trajectories, the system total energy, the radial distribution and velocity autocorrelation functions are in excellent agreement with the results from the conventional MD simulation. Our neural-network approach is generally applicable to systems with more complex interatomic interactions than the L-J potential. We demonstrated that our neural-network approach with principle of least action has the following advantages over the conventional MD method: neural network provides the entire trajectory in one fell swoop over a time period of 1 ps; timestep in the neural-network approach to action is a factor of 20 of the time step in conventional MD simulation; and even with a much larger time step, the energy conservation over the entire time domain is as good as in the conventional MD simulation. 61 We also applied neural network to explore the dataset produced from MD simulations. An RL model based on deep neural network was employed to accelerate the design of kirigami structures for a TMD material, MoS2. After training on the results from a total number of 66967 MD simulations, we observed that RL can be used to generate kirigami structures with high stretchability despite an extremely large search space consisting of millions of structures, which is otherwise impossible to determine from random search due to the sparsity of the highly stretchable case (< 5%). RL agent also reveals the importance of the location and configuration of a cut unit inside the material that will eventually lead to structures with superior mechanical properties. 7.2 Future work Going forward from our current discussion on the folding (origami) and property tuning (kirigami) techniques for 2D metamaterials, we are interested in studying the electron transport in strained origami and stretched kirigami strucutres via density functional theory calculations and quantum MD simulations. We will expect more work to explore the 2D material heterostructures in quantum realm and explain the effect of the previously-mentioned two paper arts on the electronic properties of Van der Walls 2D materials and devices. By studying the electronic behaviors such as local electron affinity, density of electronic states, band structure and ionization potential, we will have a better physical picture about the constrained and deformed origami and kirigami 2D metamaterials. As a complement to the conventional MD method, our action-based NN model will be enhanced further. Two major focuses are the applicability of our model to complex atomic systems and its computational scalability. In order to extend our NN model to scenarios other than L-J interaction, we will integrate several common force fields in our model, including the previously- mentioned TIP4P/2005 potential for condensed water, REBO and Stillinger-Weber potentials for 62 MoS2 and other 2D materials. The diffusion of water molecules and 1T to 2H phase transition of MoS2 will be studied with our network then validated using conventional MD method. In addition to the applicability of our model, we propose a time-decomposition approach to improve its scalability. In fact, our neural-network approach to principle of least action has the potential to parallelize the time axis in atomistic simulations in a way similar to the divide-and- conquer schemes used in spatial decomposition of large-scale MD simulations[141]. These schemes have made it possible to simulate multimillion-to-multibillion atom systems[142]. In the time-parallelization approach, one would split a long-time trajectory, say 1ns long, into many shorter time segments on the order of 1 ps. Using the Hamiltonian Monte Carlo approach[143] accelerated by the Boltzmann generator[144], one can sample a Markov chain of 1,000 intermediate configurations between t = 0 and t = 1 ns. Subsequently, one could use the proposed least action neural-network approach to run 1,000 simulations concurrently on many processors and generate a continuous trajectory over a time period of 1 ns. Thus, the neural-network approach would provide substantial gains in computational efficiency. 63 References [1] A. K. Geim and K. S. Novoselov, “The rise of graphene,” Nat. Mater., vol. 6, no. 3, pp. 183–191, Mar. 2007. [2] S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, and A. Kis, “2D transition metal dichalcogenides,” Nat. Rev. Mater., vol. 2, 2017. [3] T. T. Tran, K. Bray, M. J. Ford, M. Toth, and I. Aharonovich, “Quantum emission from hexagonal boron nitride monolayers,” Nat. Nanotechnol., vol. 11, no. 1, pp. 37–41, Jan. 2016. [4] M. Garcia-hernandez, J. Coleman, Z. Lin, A. Mccreary, and N. Briggs, “Fundamental study Scalable shear-exfoliation of high-quality phosphorene nano flakes with reliable electrochemical cycleability in nano batteries,” 2D Mater., vol. 3, no. 2, p. 0. [5] L. Li et al., “Direct observation of the layer-dependent electronic structure in phosphorene,” Nat. Nanotechnol., vol. 12, no. 1, pp. 21–25, Jan. 2017. [6] I. A. Ovid’Ko, “Mechanical properties of graphene,” Rev. Adv. Mater. Sci., vol. 34, pp. 1– 11, 2013. [7] M. K. Blees et al., “Graphene kirigami,” Nature, vol. 524, no. 7564, pp. 204–207, Aug. 2015. [8] C. Lee, X. Wei, J. W. Kysar, and J. Hone, “Measurement of the elastic properties and intrinsic strength of monolayer graphene,” Science., vol. 321, no. 5887, pp. 385–388, Jul. 2008. [9] M. A. Lukowski et al., “Highly active hydrogen evolution catalysis from metallic WS2 nanosheets,” Energy Environ. Sci., vol. 7, no. 8, pp. 2608–2613, Jul. 2014. [10] N. Mao, Y. Chen, D. Liu, J. Zhang, and L. Xie, “Solvatochromic effect on the photoluminescence of MoS2 monolayers,” Small, vol. 9, no. 8, pp. 1312–1315, Apr. 2013. [11] W. Han, R. K. Kawakami, M. Gmitra, and J. Fabian, “Graphene spintronics,” Nature Nanotechnology, vol. 9, no. 10. Nature Publishing Group, pp. 794–807, 01-Jan-2014. [12] H. Tang, J. Wang, H. Yin, H. Zhao, D. Wang, and Z. Tang, “Growth of polypyrrole ultrathin films on MoS2 monolayers as high-performance supercapacitor electrodes,” Adv. Mater., vol. 27, no. 6, pp. 1117–1123, 2015. [13] B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, and A. Kis, “Single-layer MoS2 transistors,” Nat. Nanotechnol., vol. 6, no. 3, pp. 147–150, Jan. 2011. [14] S. Bertolazzi, J. Brivio, and A. Kis, “Stretching and breaking of ultrathin MoS2,” ACS Nano, vol. 5, no. 12, pp. 9703–9709, Dec. 2011. [15] H. C. Diaz, R. Chaghi, Y. Ma, and M. Batzill, “Molecular beam epitaxy of the van der Waals heterostructure MoTe2 on MoS2 : phase, thermal, and chemical stability,” 2D Mater., vol. 2, no. 4, p. 044010, Nov. 2015. [16] H. J. Conley, B. Wang, J. I. Ziegler, R. F. Haglund, S. T. Pantelides, and K. I. Bolotin, “Bandgap engineering of strained monolayer and bilayer MoS2,” Nano Lett., vol. 13, no. 8, pp. 3626–3630, Aug. 2013. [17] G. H. Lee et al., “Flexible and transparent MoS2 field-effect transistors on hexagonal boron nitride-graphene heterostructures,” ACS Nano, vol. 7, no. 9, pp. 7931–7936, Sep. 2013. [18] A. Apte et al., “Structural phase transformation in strained monolayer MoWSe2 alloy,” ACS Nano, vol. 12, no. 4, pp. 3468–3476, Apr. 2018. [19] J. Shen et al., “Liquid phase exfoliation of two-dimensional materials by directly probing and matching surface tension components,” Nano Lett., vol. 15, no. 8, pp. 5449–5454, Aug. 2015. 64 [20] G. Zhou et al., “Molecular simulation of MoS2 exfoliation,” Sci. Rep., vol. 8, no. 1, p. 16761, Dec. 2018. [21] D.-G. Hwang and M. D. Bartlett, “Tunable mechanical metamaterials through hybrid kirigami structures,” Sci. Rep., vol. 8, no. 1, p. 3378, Dec. 2018. [22] Y. Tang, G. Lin, S. Yang, Y. K. Yi, R. D. Kamien, and J. Yin, “Programmable kiri-kirigami metamaterials,” Adv. Mater., vol. 29, no. 10, p. 1604262, Mar. 2017. [23] C. Lee, X. Wei, J. W. Kysar, and J. Hone, “Measurement of the elastic properties and intrinsic strength of monolayer graphene,” Science, vol. 321, no. 5887, pp. 385–388, Jul. 2008. [24] S. Zhu and T. Li, “Hydrogenation-assisted graphene origami and its application in programmable molecular mass uptake, storage, and release,” ACS Nano, vol. 8, no. 3, pp. 2864–2872, Mar. 2014. [25] Y. Cao et al., “Correlated insulator behaviour at half-filling in magic-angle graphene superlattices,” Nature, vol. 556, no. 7699, pp. 80–84, Apr. 2018. [26] J. Yin et al., “Selectively enhanced photocurrent generation in twisted bilayer graphene with van Hove singularity,” Nat. Commun., vol. 7, no. 1, pp. 1–8, Mar. 2016. [27] R. W. Havener, Y. Liang, L. Brown, L. Yang, and J. Park, “Van hove singularities and excitonic effects in the optical conductivity of twisted bilayer graphene,” Nano Lett., vol. 14, no. 6, pp. 3353–3357, Jun. 2014. [28] Y. Cao et al., “Unconventional superconductivity in magic-angle graphene superlattices,” Nature, vol. 556, no. 7699, pp. 43–50, Apr. 2018. [29] S. Carr, D. Massatt, S. Fang, P. Cazeaux, M. Luskin, and E. Kaxiras, “Twistronics: Manipulating the electronic properties of two-dimensional layered structures through their twist angle,” Phys. Rev. B, vol. 95, no. 7, p. 075420, Feb. 2017. [30] H. Zhao et al., “Fluidic flow assisted deterministic folding of van der Waals materials,” Adv. Funct. Mater., vol. 30, no. 13, p. 1908691, Mar. 2020. [31] K. Shimamura et al., “Guidelines for creating artificial neural network empirical interatomic potential from first-principles molecular dynamics data under specific conditions and its application to α-Ag2Se,” J. Chem. Phys., vol. 151, no. 12, p. 124303, Sep. 2019. [32] S. Fukushima et al., “Thermodynamic integration by neural network potentials based on first-principles dynamic calculations,” Phys. Rev. B, vol. 100, no. 21, p. 214108, Dec. 2019. [33] M. L. Piscopo, M. Spannowsky, and P. Waite, “Solving differential equations with neural networks: Applications to the calculation of cosmological phase transitions,” Phys. Rev. D, vol. 100, no. 1, p. 16002, 2019. [34] W. C. Swope, H. C. Andersen, P. H. Berens, and K. R. Wilson, “A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters,” J. Chem. Phys., vol. 76, no. 1, pp. 637–649, Jan. 1982. [35] S. Hong et al., “Computational synthesis of MoS2 layers by reactive molecular dynamics simulations: Initial sulfidation of MoO3 surfaces,” Nano Lett., vol. 17, no. 8, pp. 4866–4872, Aug. 2017. [36] V. Yamakov, D. Wolf, S. R. Phillpot, A. K. Mukherjee, and H. Gleiter, “Deformation- mechanism map for nanocrystalline metals by molecular-dynamics simulation,” Nat. Mater., vol. 3, no. 1, pp. 43–47, Dec. 2004. [37] T. A. Ho and A. Striolo, “Molecular dynamics simulation of the graphene–water interface: comparing water models,” Mol. Simul., vol. 40, no. 14, pp. 1190–1200, Nov. 2014. 65 [38] M. W. Mahoney and W. L. Jorgensen, “A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions,” J. Chem. Phys., vol. 112, no. 20, pp. 8910–8922, May 2000. [39] C. Hong, D. P. Tieleman, and Y. Wang, “Microsecond molecular dynamics simulations of lipid mixing,” Langmuir, vol. 30, no. 40, pp. 11993–12001, Oct. 2014. [40] K. Kuwata, Y. O. Kamatari, K. Akasaka, and T. L. James, “Slow conformational dynamics in the hamster prion protein † ,” Biochemistry, vol. 43, no. 15, pp. 4439–4446, Apr. 2004. [41] P. Maragakis et al., “Microsecond molecular dynamics simulation shows effect of slow loop dynamics on backbone amide order parameters of proteins,” J. Phys. Chem. B, vol. 112, no. 19, pp. 6155–6158, May 2008. [42] A. Nakano, “A space-time-ensemble parallel nudged elastic band algorithm for molecular kinetics simulation,” Comput. Phys. Commun., vol. 178, no. 4, pp. 280–289, Feb. 2008. [43] L. Bassman et al., “Active learning for accelerated design of layered materials,” npj Comput. Mater., vol. 4, no. 1, pp. 1–9, Dec. 2018. [44] S. Plimpton, “Fast parallel algorithms for short-range molecular dynamics,” J. Comput. Phys., vol. 117, no. 1, pp. 1–19, Mar. 1995. [45] V. Varshney, S. S. Patnaik, C. Muratore, A. K. Roy, A. A. Voevodin, and B. L. Farmer, “MD simulations of molybdenum disulphide (MoS2): Force-field parameterization and thermal transport behavior,” Comput. Mater. Sci., vol. 48, no. 1, pp. 101–108, 2010. [46] J. W. Jiang, H. S. Park, and T. Rabczuk, “Molecular dynamics simulations of single-layer molybdenum disulphide (MoS2): Stillinger-Weber parametrization, mechanical properties, and thermal conductivity,” J. Appl. Phys., vol. 114, no. 6, 2013. [47] J.-W. Jiang, “Parametrization of Stillinger–Weber potential based on valence force field model: Application to single-layer MoS2 and black phosphorus,” Nanotechnology, vol. 26, no. 31, p. 315706, Aug. 2015. [48] E. W. Bucholz and S. B. Sinnott, “Mechanical behavior of MoS2 nanotubes under compression, tension, and torsion from molecular dynamics simulations,” J. Appl. Phys., vol. 112, no. 12, p. 123510, Dec. 2012. [49] T. Liang, S. R. Phillpot, and S. B. Sinnott, “Parametrization of a reactive many-body potential for Mo–S systems,” Phys. Rev. B, vol. 79, no. 24, p. 245110, Jun. 2009. [50] B. Luan and R. Zhou, “Wettability and friction of water on a MoS2 nanosheet,” Appl. Phys. Lett., vol. 108, no. 13, p. 131601, Mar. 2016. [51] L. Wang, Y. Wang, J. I. Wong, T. Palacios, J. Kong, and H. Y. Yang, “Functionalized MoS2 nanosheet-based field-effect biosensor for label-free sensitive detection of cancer marker proteins in solution,” Small, vol. 10, no. 6, pp. 1101–1105, Mar. 2014. [52] A. Smolyanitsky, B. I. Yakobson, T. A. Wassenaar, E. Paulechka, and K. Kroenlein, “A MoS2-based capacitive displacement sensor for DNA sequencing,” ACS Nano, vol. 10, no. 9, pp. 9009–9016, Sep. 2016. [53] M. C. Gordillo, G. Nagy, and J. Martí, “Structure of water nanoconfined between hydrophobic surfaces,” J. Chem. Phys., vol. 123, no. 5, p. 054707, Aug. 2005. [54] I. C. Bourg and C. I. Steefel, “Molecular dynamics simulations of water structure and diffusion in silica nanopores,” J. Phys. Chem. C, vol. 116, no. 21, pp. 11556–11564, May 2012. [55] Q. Gao et al., “Effect of adsorbed alcohol layers on the behavior of water molecules confined in a graphene nanoslit: A molecular dynamics study,” Langmuir, vol. 33, no. 42, pp. 11467–11474, Oct. 2017. 66 [56] P. Kim, T. S. Wong, J. Alvarenga, M. J. Kreder, W. E. Adorno-Martinez, and J. Aizenberg, “Liquid-infused nanostructured surfaces with extreme anti-ice and anti-frost performance,” ACS Nano, vol. 6, no. 8, pp. 6569–6577, 2012. [57] J. S. Palmer, S. Sivaramakrishnan, P. S. Waggoner, and J. H. Weaver, “Particle aggregation on dewetting solid water films,” Surf. Sci., vol. 602, no. 13, pp. 2278–2283, 2008. [58] Y. Zheng et al., “Directional water collection on wetted spider silk,” Nature, vol. 463, no. 7281, pp. 640–643, 2010. [59] T. R. Jensen et al., “Water in contact with extended hydrophobic surfaces: Direct evidence of weak dewetting,” Phys. Rev. Lett., vol. 90, no. 8, p. 4, 2003. [60] A. Kayal and A. Chandra, “Wetting and dewetting of narrow hydrophobic channels by orthogonal electric fields: Structure, free energy, and dynamics for different water models,” J. Chem. Phys., vol. 143, no. 22, p. 224708, Dec. 2015. [61] J. Zhang et al., “Molecular insight into nanoscale water films dewetting on modified silica surfaces,” Phys. Chem. Chem. Phys., vol. 17, no. 1, pp. 451–458, Jan. 2015. [62] C.-J. Shih, S. Lin, M. S. Strano, and D. Blankschtein, “Understanding the stabilization of liquid-phase-exfoliated graphene in polar solvents: Molecular dynamics simulations and kinetic theory of colloid aggregation,” J. Am. Chem. Soc., vol. 132, no. 41, pp. 14638–14648, Oct. 2010. [63] J.-M. Yun et al., “Exfoliated and partially oxidized MoS2 nanosheets by one-pot reaction for efficient and stable organic solar cells,” Small, vol. 10, no. 12, pp. 2319–2324, 2014. [64] M. Wang et al., “Surface tension components ratio: An efficient parameter for direct liquid phase exfoliation,” ACS Appl. Mater. Interfaces, vol. 9, no. 10, pp. 9168–9175, 2017. [65] D. Hanlon et al., “Production of molybdenum trioxide nanosheets by liquid exfoliation and their application in high-performance supercapacitors,” Chem. Mater., vol. 26, no. 4, pp. 1751–1763, 2014. [66] H. Chen et al., “Atomically precise, custom-design origami graphene nanostructures,” Science., vol. 365, no. 6457, pp. 1036–1040, Sep. 2019. [67] S. Cranford, D. Sen, and M. J. Buehler, “Meso-origami: Folding multilayer graphene sheets,” Appl. Phys. Lett., vol. 95, no. 12, p. 123121, Sep. 2009. [68] J. Mu, C. Hou, H. Wang, Y. Li, Q. Zhang, and M. Zhu, “Origami-inspired active graphene- based paper for programmable instant self-folding walking devices,” Sci. Adv., vol. 1, no. 10, p. e1500533, Nov. 2015. [69] P. Y. Huang et al., “Grains and grain boundaries in single-layer graphene atomic patchwork quilts,” Nature, vol. 469, no. 7330, pp. 389–392, Jan. 2011. [70] P. Johari and V. B. Shenoy, “Tuning the electronic properties of semiconducting transition metal dichalcogenides by applying mechanical strains,” ACS Nano, vol. 6, no. 6, pp. 5449– 5456, Jun. 2012. [71] W. Wu et al., “Piezoelectricity of single-atomic-layer MoS2 for energy conversion and piezotronics,” Nature, vol. 514, no. 7523, pp. 470–474, Oct. 2014. [72] Y. Yang et al., “Brittle Fracture of 2D MoSe2,” Adv. Mater., vol. 29, no. 2, p. 1604201, Jan. 2017. [73] L. Cai et al., “Chemically derived kirigami of WSe2,” J. Am. Chem. Soc., vol. 140, no. 35, pp. 10980–10987, Sep. 2018. [74] W. Zheng et al., “Kirigami-inspired highly stretchable nanoscale devices using multidimensional deformation of monolayer MoS2,” Chem. Mater., vol. 30, no. 17, pp. 6063–6070, Sep. 2018. [75] A. Föppl, Vorlesungen über technische Mechanik. Leipzig, B. G. Teubner, 1907. [76] T. von Kármán, Festigkeitsproblem im Maschinenbau. Encyk. D. Math. Wiss., 1910. 67 [77] T. Liang, S. R. Phillpot, and S. B. Sinnott, “Erratum: Parametrization of a reactive many- body potential for Mo-S systems,” Phys. Rev. B - Condens. Matter Mater. Phys., vol. 85, no. 19, pp. 2–3, 2012. [78] K. Liu et al., “Elastic properties of chemical-vapor-deposited monolayer MoS2, WS2, and their bilayer Heterostructures,” Nano Lett., vol. 14, no. 9, pp. 5097–5103, Sep. 2014. [79] A. B. Churnside et al., “Routine and timely sub-pico Newton force stability and precision for biological applications of atomic force microscopy,” Nano Lett., vol. 12, no. 7, pp. 3557–3561, Jul. 2012. [80] M. Huang, T. A. Pascal, H. Kim, W. A. Goddard, and J. R. Greer, “Electronic−mechanical coupling in graphene from in situ nanoindentation experiments and multiscale atomistic simulations,” Nano Lett., vol. 11, no. 3, pp. 1241–1246, Mar. 2011. [81] J. A. Stewart and D. E. Spearot, “Atomistic simulations of nanoindentation on the basal plane of crystalline molybdenum disulfide (MoS2),” Model. Simul. Mater. Sci. Eng., vol. 21, no. 4, 2013. [82] X. Tan, J. Wu, K. Zhang, X. Peng, L. Sun, and J. Zhong, “Nanoindentation models and Young’s modulus of monolayer graphene: A molecular dynamics study,” Appl. Phys. Lett., vol. 102, no. 7, 2013. [83] T. C. Shyu et al., “A kirigami approach to engineering elasticity in nanocomposites through patterned defects,” Nat. Mater., vol. 14, no. 8, pp. 785–789, Aug. 2015. [84] P. Z. Hanakata, Z. Qi, D. K. Campbell, and H. S. Park, “Highly stretchable MoS2 kirigami,” Nanoscale, vol. 8, no. 1, pp. 458–463, 2016. [85] W. S. McCulloch and W. Pitts, “A logical calculus of the ideas immanent in nervous activity,” Bull. Math. Biophys., vol. 5, no. 4, pp. 115–133, Dec. 1943. [86] F. Rosenblatt, “The perceptron: A probabilistic model for information storage and organization in the brain,” Psychol. Rev., vol. 65, no. 6, pp. 386–408, 1958. [87] J. J. Hopfield, “Neural networks and physical systems with emergent collective computational abilities,” Proc. Natl. Acad. Sci., vol. 79, no. 8, pp. 2554–2558, Apr. 1982. [88] J. Behler, “Perspective: Machine learning potentials for atomistic simulations,” J. Chem. Phys., vol. 145, no. 17, p. 170901, Nov. 2016. [89] J. Lim, S. Ryu, J. W. Kim, and W. Y. Kim, “Molecular generative model based on conditional variational autoencoder for de novo molecular design,” J. Cheminform., vol. 10, no. 1, p. 31, Dec. 2018. [90] P. Rajak, A. Krishnamoorthy, A. Nakano, P. Vashishta, and R. Kalia, “Structural phase transitions in a MoWSe2 monolayer: Molecular dynamics simulations and variational autoencoder analysis,” Phys. Rev. B, vol. 100, no. 1, p. 014108, Jul. 2019. [91] J. Liu et al., “Boltzmann machine modeling of layered MoS2 synthesis on a quantum annealer,” Comput. Mater. Sci., vol. 173, p. 109429, Feb. 2020. [92] D. Voiry, J. Yang, and M. Chhowalla, “Recent strategies for improving the catalytic activity of 2D TMD nanosheets toward the hydrogen evolution reaction,” Adv. Mater., vol. 28, no. 29, pp. 6197–6206, Aug. 2016. [93] L. Xie and X. Cui, “Manipulating spin-polarized photocurrents in 2D transition metal dichalcogenides,” Proc. Natl. Acad. Sci., vol. 113, no. 14, pp. 3746–3750, Apr. 2016. [94] E. Okogbue et al., “Multifunctional two-dimensional PtSe2-layer kirigami conductors with 2000% stretchability and metallic-to-semiconducting tunability,” Nano Lett., vol. 19, no. 11, pp. 7598–7607, Nov. 2019. [95] R. M. Neville, F. Scarpa, and A. Pirrera, “Shape morphing kirigami mechanical metamaterials,” Sci. Rep., vol. 6, no. 1, p. 31067, Nov. 2016. 68 [96] Z. G. Nicolaou and A. E. Motter, “Mechanical metamaterials with negative compressibility transitions,” Nat. Mater., vol. 11, no. 7, pp. 608–613, Jul. 2012. [97] J. N. Grima, D. Attard, and R. Gatt, “Truss-type systems exhibiting negative compressibility,” Phys. status solidi, vol. 245, no. 11, pp. 2405–2414, Nov. 2008. [98] P. Z. Hanakata, E. D. Cubuk, D. K. Campbell, and H. S. Park, “Accelerated search and design of stretchable graphene kirigami using machine learning,” Phys. Rev. Lett., vol. 121, no. 25, p. 255304, Dec. 2018. [99] S. Zhu, Y. Huang, and T. Li, “Extremely compliant and highly stretchable patterned graphene,” Appl. Phys. Lett., vol. 104, no. 17, p. 173103, Apr. 2014. [100] Z. Liu, D. Zhu, S. P. Rodrigues, K.-T. Lee, and W. Cai, “Generative model for the inverse design of metasurfaces,” Nano Lett., vol. 18, no. 10, pp. 6570–6576, Oct. 2018. [101] Z. Zhou, X. Li, and R. N. Zare, “Optimizing chemical reactions with deep reinforcement learning,” ACS Cent. Sci., vol. 3, no. 12, pp. 1337–1344, Dec. 2017. [102] Z. Zhou, S. Kearnes, L. Li, R. N. Zare, and P. Riley, “Optimization of molecules via deep reinforcement learning,” Sci. Rep., vol. 9, no. 1, pp. 1–10, Dec. 2019. [103] J. L. F. Abascal and C. Vega, “A general purpose model for the condensed phases of water: TIP4P/2005,” J. Chem. Phys., vol. 123, no. 23, p. 234505, Dec. 2005. [104] W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, “Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids,” J. Am. Chem. Soc., vol. 118, no. 45, pp. 11225–11236, Jan. 1996. [105] V. Sresht, A. Govind Rajan, E. Bordes, M. S. Strano, A. A. H. Pádua, and D. Blankschtein, “Quantitative modeling of MoS2–solvent interfaces: predicting contact angles and exfoliation performance using molecular dynamics,” J. Phys. Chem. C, vol. 121, no. 16, pp. 9022–9031, Apr. 2017. [106] A. P. S. Gaur, S. Sahoo, M. Ahmadi, S. P. Dash, M. J. F. Guinel, and R. S. Katiyar, “Surface energy engineering for tunable wettability through controlled synthesis of MoS2,” Nano Lett., vol. 14, no. 8, pp. 4314–4321, Aug. 2014. [107] H. Peng, Z.-H. Yang, J. P. Perdew, and J. Sun, “Versatile van der Waals density functional based on a meta-generalized gradient approximation,” Phys. Rev. X, vol. 6, no. 4, p. 041005, Oct. 2016. [108] P.-G. de. Gennes, F. Brochard-Wyart, and D. Quéré, Capillarity and wetting phenomena : drops, bubbles, pearls, waves, Illustrate. Springer Science & Business Media, 2013. [109] S. J. Stuart, A. B. Tutein, and J. A. Harrison, “A reactive potential for hydrocarbons with intermolecular interactions,” J. Chem. Phys., 2000. [110] M. Ma, G. Tocci, A. Michaelides, and G. Aeppli, “Fast diffusion of water nanodroplets on graphene,” Nat. Mater., vol. 15, no. 1, pp. 66–71, Jan. 2016. [111] C. Jarzynski, “Nonequilibrium equality for free energy differences,” Phys. Rev. Lett., vol. 78, no. 14, pp. 2690–2693, Apr. 1997. [112] G. F. Schneider et al., “Tailoring the hydrophobicity of graphene for its use as nanopores for DNA translocation,” Nat. Commun., vol. 4, no. 1, pp. 1–7, Oct. 2013. [113] B. Liu et al., “Thermal transport in a graphene-MoS2 bilayer heterostructure: A molecular dynamics study,” RSC Adv., vol. 5, no. 37, pp. 29193–29200, 2015. [114] F. H. Stillinger and T. A. Weber, “Computer simulation of local order in condensed phases of silicon,” Phys. Rev. B, vol. 31, no. 8, pp. 5262–5271, Apr. 1985. [115] L. D. Landau and E. M. Lifshitz, Course of Theoretical Physics Volume 1: Mechanics, 3rd ed. Oxford: Butterworth-Heinemann, 1982. [116] L. D. Landau and E. M. Lifshitz, Course of Theoretical Physics Volume 2: The Classical Theory of Fields, 4th ed. Oxford: Butterworth-Heinemann, 1987. 69 [117] R. P. Feynman, “The principle of least action in quantum mechanics,” in Feynman’s Thesis — A New Approach to Quantum Theory, World Scientific, 2005, pp. 1–69. [118] V. García-Morales, J. Pellicer, and J. A. Manzanares, “Thermodynamics based on the principle of least abbreviated action: Entropy production in a network of coupled oscillators,” Ann. Phys. (N. Y)., vol. 323, no. 8, pp. 1844–1858, Aug. 2008. [119] R. P. Feynman, R. B. Leighton, M. Sands, and E. M. Hafner, “The Feynman Lectures on Physics; Vol. I,” Am. J. Phys., vol. 33, no. 9, pp. 750–752, Sep. 1965. [120] H. Goldstein, Classical Mechanics, 3rd ed. Boston: Addison-Wesley, 2002. [121] A. Rahman, “Correlations in the motion of atoms in liquid argon,” Phys. Rev., vol. 136, no. 2A, pp. A405–A411, Oct. 1964. [122] A. C. Rajan et al., “Machine-learning-assisted accurate band gap predictions of functionalized mxene,” Chem. Mater., vol. 30, no. 12, pp. 4031–4038, Feb. 2018. [123] J. Wang, X. Yang, Z. Zeng, X. Zhang, X. Zhao, and Z. Wang, “New methods for prediction of elastic constants based on density functional theory combined with machine learning,” Comput. Mater. Sci., vol. 138, pp. 135–148, Oct. 2017. [124] Y. Iwasaki et al., “Machine-learning guided discovery of a new thermoelectric material,” Sci. Rep., vol. 9, no. 1, pp. 1–7, Dec. 2019. [125] T. Wang, C. Zhang, H. Snoussi, and G. Zhang, “Machine learning approaches for thermoelectric materials research,” Adv. Funct. Mater., vol. 30, no. 5, p. 1906041, Jan. 2020. [126] L. Onsager and S. MacHlup, “Fluctuations and irreversible processes,” Phys. Rev., vol. 91, no. 6, pp. 1505–1512, Sep. 1953. [127] L. Onsager, “Reciprocal relations in rrreversible processes. I.,” Phys. Rev., vol. 37, no. 4, pp. 405–426, Feb. 1931. [128] L. Onsager, “Reciprocal relations in irreversible processes. II.,” Phys. Rev., vol. 38, no. 12, pp. 2265–2279, Dec. 1931. [129] D. E. Shaw et al., “Atomic-level characterization of the structural dynamics of proteins,” Science., vol. 330, no. 6002, pp. 341–346, Oct. 2010. [130] C. Krzeminski, Q. Brulin, V. Cuny, E. Lecat, E. Lampin, and F. Cleri, “Molecular dynamics simulation of the recrystallization of amorphous Si layers: Comprehensive study of the dependence of the recrystallization velocity on the interatomic potential,” J. Appl. Phys., vol. 101, no. 12, p. 123506, Jun. 2007. [131] M. P. Allen and D. J. Tildesley, Computer simulation of liquids: Second edition, 2nd ed. Oxford: Oxford University Press, 2017. [132] K. Yasue, “The role of the Onsager–Machlup Lagrangian in the theory of stationary diffusion process,” J. Math. Phys., vol. 20, no. 9, pp. 1861–1864, Sep. 1979. [133] D. Passerone and M. Parrinello, “Action-derived molecular dynamics in the study of rare events,” Phys. Rev. Lett., vol. 87, no. 10, pp. 108302/1-108302/4, Sep. 2001. [134] I. H. Lee, J. Lee, and S. Lee, “Kinetic energy control in action-derived molecular dynamics simulations,” Phys. Rev. B - Condens. Matter Mater. Phys., vol. 68, no. 6, p. 064303, Aug. 2003. [135] J. Lee, I. H. Lee, I. S. Joung, J. Lee, and B. R. Brooks, “Finding multiple reaction pathways via global optimization of action,” Nat. Commun., vol. 8, no. 1, pp. 1–8, May 2017. [136] T. Dozat, “Incorporating Nesterov Momentum into Adam,” ICLR Workshop, no. 1. pp. 2013–2016, 2016. [137] R. Frostig, M. J. Johnson, and C. Leary, “Compiling machine learning programs via high- level tracing,” SysML. 2018. [138] J. G. Kirkwood and E. M. Boggs, “The radial distribution function in liquids,” J. Chem. Phys., vol. 10, no. 6, pp. 394–402, Jun. 1942. 70 [139] B. J. Alder and T. E. Wainwright, “Decay of the velocity autocorrelation function,” Phys. Rev. A, vol. 1, no. 1, pp. 18–21, Jan. 1970. [140] D. J. A. Hills, A. M. Grütter, and J. J. Hudson, “An algorithm for discovering Lagrangians automatically from data,” PeerJ Comput. Sci., vol. 1, no. 11, p. e31, Nov. 2015. [141] F. Shimojo, R. K. Kalia, A. Nakano, and P. Vashishta, “Embedded divide-and-conquer algorithm on hierarchical real-space grids: parallel molecular dynamics simulation based on linear-scaling density functional theory,” Comput. Phys. Commun., vol. 167, no. 3, pp. 151– 164, May 2005. [142] A. Nakano et al., “A divide-and-conquer/cellular-decomposition framework for million-to- billion atom simulations of chemical reactions,” Comput. Mater. Sci., vol. 38, no. 4, pp. 642–652, Feb. 2007. [143] K. D. Dang, M. Quiroz, R. Kohn, M. N. Tran, and M. Villani, “Hamiltonian monte carlo with energy conserving subsampling,” 2019. [144] F. Noé, S. Olsson, J. Köhler, and H. Wu, “Boltzmann generators: sampling equilibrium states of many-body systems with deep learning,” Science., vol. 365, no. 6457, Sep. 2019.

Abstract (if available)

Abstract In recent years, two-dimensional (2D) materials have become the spotlight of nanoscience owing to their exotic properties. Massive experimental and theoretical research has been triggered to understand the synthesis of 2D materials via exfoliation and to explore the controllability of their properties via structure engineering. As a synthesis technique, liquid-phase exfoliation is of great significance due to its simplicity and high yield. Researchers have devoted many efforts to explain and improve the synthesis of 2D materials. However, understanding and optimizing the synthesis experiment at nanoscale remains challenging. Structure engineering methods such as origami and kirigami offer efficient approaches to designing and tuning novel 2D metamaterial structures and properties. Folding and stacking of atomically-thin nanosheets can lead to intriguing new physical phenomena such as bandgap tuning, metal-insulator transition, Van Hove singularity and superconductivity. Nevertheless, achieving scalable folding and reasonable designing of 2D materials with high spatial precision and desired properties have been difficult tasks. To better understand the synthesis of 2D materials via liquid-phase exfoliation, to develop a general method for origami and kirigam structure designing, we first perform molecular dynamics simulations to study the interlayer coupling and delamination of a 2D material in the presence of liquid solvents and then propose a scalable deterministic technique to fold the synthesized 2D metamaterials with a pre-defined configuration. We investigate the mechanical properties of kirigami structures of a 2D material (i.e. MoS₂) and discuss the underlying mechanism and energetics of the process. Finally, we demonstrate how neural-network based models can accelerate the design of kirigami 2D materials and address foundational problems in molecular dynamics.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Large-scale molecular dynamics simulations of nano-structured materials

PDF

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

PDF

High throughput computational framework for synthesis and accelerated discovery of dielectric polymer materials using polarizable reactive molecular dynamics and graph neural networks

PDF

Physics informed neural networks and electrostrictive cavitation in water

PDF

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

PDF

Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations

PDF

Simulation and machine learning at exascale

PDF

Shock-induced poration, cholesterol flip-flop and small interfering RNA transfection in a phospholipid membrane: multimillion atom, microsecond molecular dynamics simulations

PDF

Molecular dynamics simulations of nanostructures

PDF

Thermal properties of silicon carbide and combustion mechanisms of aluminum nanoparticle

PDF

Light Emission from Carbon Nanotubes and Two-Dimensional Materials

PDF

Shock sensitivity of energetic material and nanometric damage mechanisms in silica glass

PDF

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

PDF

Coulomb interactions and superconductivity in low dimensional materials

PDF

Dynamics of wing cracks and nanoscale damage in silica glass

PDF

Reactive and quantum molecular dynamics study of materials: from oxidation and amorphization to interfacial charge transfer

PDF

Theoretical modeling of nanoscale systems: applications and method development

PDF

2D layered materials: fundamental properties and device applications

$Multimillion-to-billion atom molecular dynamics simulations of deformation, damage, nanoindentation, and fracture in silica glass and energetic materials$

PDF

Multimillion-to-billion atom molecular dynamics simulations of deformation, damage, nanoindentation, and fracture in silica glass and energetic materials

PDF

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

Asset Metadata

Creator Wang, Beibei (author)

Core Title Neural network for molecular dynamics simulation and design of 2D materials

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Physics

Publication Date 07/10/2020

Defense Date 04/17/2020

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag 2D materials,molecular dynamics,neural network,OAI-PMH Harvest

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Kalia, Rajiv (committee chair), Haas, Stephan (committee member), Kumar, Satish (committee member), Nakano, Aiichiro (committee member), Vashishta, Priya (committee member)

Creator Email beibeiw@usc.edu,wbb.ustc@gmail.com

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

Unique identifier UC11663910

Identifier etd-WangBeibei-8646.pdf (filename),usctheses-c89-327095 (legacy record id)

Legacy Identifier etd-WangBeibei-8646.pdf

Dmrecord 327095

Document Type Dissertation

Rights Wang, Beibei

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA