ReaxFF Potential Research Articles

Molecular Dynamics Simulations of Binary Phosphate Glass Using the ReaxFF Potential.

Due to the importance of the understanding of dissolution behavior of phosphate-based bioglasses (PBGs) in different biomedical applications, binary sodium and calcium phosphate glasses have been simulated for the first time using a newly developed ReaxFF force field and a standard melt-quench method with the LAMMPS classical molecular dynamics software. The partial radial distribution function of P-O within the first coordination shell indicated two distinct peaks corresponding to phosphorus bonding to NBO and BO, respectively, at distances consistent with those observed experimentally and a P-O coordination number of 4.0. Angular distribution functions were consistent with the experimental data. The calculated network connectivities are in good agreement with experimental data, and the detailed Qn distributions are broader than would be expected.

Molecular Dynamics of Catalyst-Free Edge Elongation of Boron Nitride Nanotubes Coaxially Grown on Single-Walled Carbon Nanotubes.

Recent advances in low-dimensional materials have enabled the synthesis of single-walled carbon nanotubes encapsulated in hexagonal boron nitride (BN) nanotubes (SWCNT@BNNT), creating one-dimensional van der Waals (vdW) heterostructures. However, controlling the quality and crystallinity of BNNT on the surface of SWCNTs using chemical vapor deposition (CVD) remains a challenge. To better understand the growth mechanism of the BNNT in SWCNT@BNNT, we conducted molecular dynamics (MD) simulations using empirical potentials. The simulation results suggest that spontaneous BN nucleation is unlikely to occur on the outer surface of the SWCNT when we assume only vdW interaction between the BN and SWCNT layers. However, we observe the elongation of the BNNT when a short BNNT is provided as a seed nucleus on the SWCNT. This grown BNNT structure, with its sharply cut edges, aligns with experimental observations made using transmission electron microscopy (TEM). Moreover, the edge-reconstruction process favors zigzag B edges, which exhibit low edge energy according to the ReaxFF potential. Our simulation successfully provides insights into the catalyst-free growth process of this one-dimensional vdW heterostructure.

Modification of a ReaxFF potential at short range for energetic materials

The ReaxFF can describe the properties of energetic materials (EMs) at equilibrium state, but does not work properly in simulating high-energy particle irradiation process because of its weak short-range interaction. In this paper, a modification was made for such a potential by connecting Ziegler-Biersack-Littmark (ZBL) potential to ReaxFF-lg through comparing to Density Functional Theory (DFT) results to accurately describe short-range interactions. After modification, the newly fitted ReaxFF-lg/ZBL potential predicts better the equation of state for EMs In displacement cascade simulations, comparing to results from ab initio molecular dynamics (AIMD), ReaxFF-lg/ZBL presented the similar transferred energy from a primary knock-on atom to surrounding atoms, better than the original ReaxFF-lg potential. Further large-scale displacement cascade simulations indicated ReaxFF-lg/ZBL could be applied for cascade simulations with PKA energy from less than 1 keV to high energy (e.g. 35 keV) cases, which is suitable for effectively simulating high-energy displacement cascades in EMs using molecular dynamics method.

Physical, chemical and morphological evolution of incipient soot obtained from molecular dynamics simulation of acetylene pyrolysis

Incipient soot particles obtained from a series of reactive molecular dynamics simulations were studied to understand the evolution of physical, chemical, and morphological properties of incipient soot. Reactive molecular dynamics simulations of acetylene pyrolysis were performed using ReaxFF potential at 1350, 1500, 1650, and 1800 K. A total of 3324 incipient soot particles were extracted from the simulations at various stages of development. Features such as the number of carbon and hydrogen atoms, number of ring structures, mass, C/H ratio, radius of gyration, surface area, volume, atomic fractal dimension, and density were calculated for each particle. The calculated values of density and C/H ratio matched well with experimental values reported in the literature. Based on the calculated features, the particles were classified in two types: type 1 and type 2 particles. It was found that type 1 particles show significant morphological evolution while type 2 particles undergo chemical restructuring without any significant morphological change. The particle volume was found to be well-correlated with the number of carbon atoms in both type 1 and type 2 particles, whereas surface area was found to be correlated with the number of carbon atoms only for type 1 particles. A correlation matrix comparing the level of correlation between any two features for both type 1 and type 2 particles was created. Finally, based on the calculated statistics, a set of correlations among various physical and morphological parameters of incipient soot was proposed.

Internal Structure of Incipient Soot from Acetylene Pyrolysis Obtained via Molecular Dynamics Simulations.

A series of reactive molecular dynamics simulations is used to study the internal structure of incipient soot particles obtained from acetylene pyrolysis. The simulations were performed using the ReaxFF potential at four different temperatures. The resulting soot particles are cataloged and analyzed to obtain statistics of their mass, volume, density, C/H ratio, number of cyclic structures, and other features. A total of 3324 incipient soot particles were analyzed in this study. Based on their structural characteristics, the incipient soot particles are classified into two classes, termed type 1 and type 2 incipient soot particles in this work. The radial distribution of density, cyclic (5-, 6-, or 7-member rings) structures, and C/H ratio inside the particles revealed a clear difference in the internal structure between type 1 and type 2 particles. These classes were further found to be well represented by the size of the particles, with smaller particles in type 1 and larger particles in type 2. The radial distributions of ring structures, density, and the C/H ratio indicated the presence of a dense core region in type 2 particles. In contrast, no clear evidence of the presence of a core was found in type 1 particles. In type 2 incipient soot particles, the boundary between the core and shell was found to be around 50-60% of the particle's radius of gyration.

Morphological, structural and electronic properties of zinc oxide nano structures deposited on ZnO substrate layer – A theoretical perspective

The work investigates the deposition of zinc and oxygen atoms on a Zinc terminated ZnO substrate, maintaining their stoichiometric ratio (O:Zn) at 0.8, under varying incident energies (0.5 eV and 5 eV) and temperatures (300 K and 700 K) using the ReaxFF potential in Molecular Dynamics (MD) simulation. The study explores the impact of deposition conditions on growth mechanisms, film morphology, and structural parameters such as Radial Distribution Function (RDF) and Coordination Number (CN). RDF analysis highlights disordered atomic arrangements and diverse bonding configurations within the film, with CN showing deviations from perfect crystalline ZnO, indicating non-uniform, non-crystalline, and nano behaviour. Post-heating the films to temperatures approaching the ZnO melting point results in enhanced crystallinity, as evidenced by distinct, dominant, and multiple RDF peaks. Density Functional Theory (DFT) approach of the deposited structures elucidates the formation of ZnO nanoclusters-(ZnO)3 and (ZnO)4 and investigates their electronic structure properties. DFT analysis of the post-heated films reveals the formation of ZnO monomers (ZnO), indicating an approaching crystalline behaviour, yet the films exhibit non-epitaxial and non-uniform characteristics. DFT calculations were instrumental in validating the observed structural behaviours in the deposited films. The study underscores the significance of stoichiometric ratio in film deposition and demonstrates the temperature-dependent evolution of structural properties towards crystallinity, shedding light on the intricate dynamics of ZnO film growth offering tailored control over film properties for enhanced performance in practical applications such as optoelectronics and gas sensing.

Atomistic insights into chemical vapor deposition process of preparing silicon carbide materials using ReaxFF-MD simulation

SiC is a kind of key advanced coating material used widely. The in-depth understanding of the CVD process of preparing SiC materials at the atomic scale is beneficial for studying the CVD mechanism and then optimizing the materials. Here, a CVD model of preparing SiC materials is established based on ReaxFF-MD considering the whole pyrolysis-adsorption-desorption-growth mechanism. The accuracy of the ReaxFF potential and its parameters is verified by DFT at first. Then an energy-based Monte Carlo (EMC)-MD-CVD model is proposed based on the bonding formation energy originating from DFT. It is used to simulate the CVD process of preparing SiC materials. The CVD of MTS is studied by deposition efficiency and deposited film components at different temperatures, and it can be found it in good agreement with experimental results. Finally, the crystallization of SiC film prepared by the CVD of MTS is studied by energy evolution and atomic structure evolution. The results are helpful to understand the CVD mechanism for preparing coatings and it can be used to optimize the CVD parameters of preparing SiC materials in the future.

Degradation behavior of DLC thin film under gamma-rays irradiation using molecular dynamics simulation

γ-rays have high energy, thus, when irradiated on the materials, lead to the degradation process of the materials which means the durability and efficiency of the system decrease. Therefore, analysis of the effect of γ-rays on the materials is important, despite this, little work has so far been undertaken to predict the irradiation effect. Furthermore, because the previous theoretical equations have many assumptions, the accuracy of the prediction is low. For these reasons, molecular dynamics simulation was used to analyze and predict the effect of γ-rays on the DLC thin film in this study. This study investigates the degradation process, bond distribution, and radical creation during γ-rays irradiated on the DLC thin film. The primary knock-on atoms (PKA) method was applied which can be mimicked first excited atoms by γ-rays. Also, to analyze the chemical change of the irradiated DLC, ReaxFF potential was selected. The dose rate condition covered various environments where the γ-rays were irradiated such as space, nuclear, medical, etc. The degradation mechanisms of DLC were different by the dose rate. For low dose rate conditions, the energy of excited atoms was converted by physical collisions between atoms. However, in the case of high dose rate conditions, radical production rapidly increased, resulting in physical and chemical energy dissipation. Within this paper, we present a fundamental numerical analysis of the degradation process that combines the PKA method and MD simulation to achieve improved accuracy. These findings could provide fundamental insight into the degradation phenomenon of surface-protective coatings for various applications with γ-ray irradiation environment.

Characterization of iron(III) in aqueous and alkaline environments with abinitio and ReaxFF potentials.

The iron(III) complexes [Fe(H2O)n(OH)m]3-m (n + m = 5, 6, m ≤ 3) and corresponding proton transfer reactions are studied with total energy calculations, the nudged elastic band (NEB) method, and molecular dynamics (MD) simulations using abinitio and a modification of reactive force field potentials, the ReaxFF-AQ potentials, based on the implementation according to Böhm et al. [J. Phys. Chem. C 120, 10849-10856 (2016)]. Applying abinitio potentials, the energies for the reactions [Fe(H2O)n(OH)m]3-m + H2O → [Fe(H2O)n-1(OH)m+1]2-m + H3O+ in a gaseous environment are in good agreement with comparable theoretical results. In an aqueous (aq) or alkaline environment, with the aid of NEB computations, respective minimum energy paths with energy barriers of up to 14.6 kcal/mol and a collective transfer of protons are modeled. Within MD simulations at room temperature, a permanent transfer of protons around the iron(III) ion is observed. The information gained concerning the geometrical and energetic properties of water and the [Fe(H2O)n(OH)m]3-m complexes from the abinitio computations has been used as reference data to optimize parameters for the O-H-Fe interaction within the ReaxFF-AQ approach. For the optimized ReaxFF-AQ parameter set, the statistical properties of the basic water model, such as the radial distribution functions and the proton hopping functions, are evaluated. For the [Fe(H2O)n(OH)m]3-m complexes, it was found that while geometrical and energetic properties are in good agreement with the abinitio data for gaseous environment, the statistical properties as obtained from the MD simulations are only partly in accordance with the abinitio results for the iron(III) complexes in aqueous or alkaline environments.

Molecular mechanism of cleavage of SARS-CoV-2 spike protein by plasma generated RONS

Recently, it is been shown that cold atmospheric pressure plasmas Cold Atmospheric Plasma effectively inactivate the 2019-nCoV virus. Despite this promising finding, the precise mechanism of this inactivation remains unclear due to the limited number of studies conducted on the subject. Consequently, this paper focuses on the spike protein, a crucial part of the novel coronavirus, and the various reactive oxygen and nitrogen species (RONS) generated by the plasma. The study employs reactive molecular dynamics simulation and ReaxFF potential to explore the reactions between the spike protein molecules and different reactive oxygen nitrogen species (including H2O2, OH, O, O3, HOONO, and 1O2). The findings suggest that when a single RONS interacts with the spike protein, 1O2 and HOONO have the most potent ability to sever the spike protein. Additionally, the combined effect of long-lived and short-lived RONS presents a more potent decomposition impact.

Deformation behaviors of hydrogen filled boron nitride and boron nitride - carbon nanotubes: Molecular dynamics simulations of proposed materials for hydrogen storage, gas sensing, and radiation shielding

Boron nitride nanotubes (BNNTs) have been extensively studied for hydrogen storage, sensing, and radiation shielding applications. The deformation behaviors of pristine, defective, and H2-encapsulated closed-end armchair BNNTs and BNNT-carbon nanotube (CNT) heterostructures are investigated using molecular dynamics employing the ReaxFF interatomic potential and a combination of Tersoff, adaptive intermolecular reactive empirical bond order, and Lennard-Jones potentials. The calculated elastic properties for pristine BNNTs and CNTs are in excellent agreement with published reliable data. The nanotube length and diameter, defects, entrapped hydrogen content, and temperature all affect the mechanical properties significantly. The ReaxFF potential predicts H2 molecule adsorption on B atoms of undeformed BNNTs whilst under strain, the H2 molecules show preference for the N atoms. Hydrogen adsorption leads to lowering of the elastic and fracture strength of BNNTs. The ReaxFF potential is shown to be a reliable forcefield for investigating the interaction of H2 molecules with BNNTs and BNNT-CNT heterotubes.

Interactions of graphene with oxidants in a mixed atmosphere: synergistic effects of O2/H2O and O2/CO2 on gasification reactivity and kinetics.

The gasification of carbon with O2, CO2, and H2O oxidants plays an important role in several energy-based applications. As most of the industrial gasification processes are conducted under mixed-atmosphere conditions, the oxidation of carbon in binary oxidant mixtures becomes crucially important. Using reactive force-field (ReaxFF) potentials, extensive MD simulations were carried out on the oxidation behavior of graphene in mixed O2/H2O and O2/CO2 environments for a range of gas compositions and temperatures. A graphene sheet with a line defect comprising of eight and four-membered rings was used as the starting carbon structure. In addition to enhanced carbon gasification with oxygen additions, MD simulations showed synergistic interactions between different oxidants and their net influence on the overall reactivities. The gasification levels achieved under the binary system were higher than the linear combination of contributions from individual oxidants. The addition of ∼40% O2 in the binary mix was identified as the region with the highest reactivity during the initial stages of gasification. The oxidation reactions with oxygen were found to start instantaneously in the presence of H2O or CO2 instead of the usual initial delay. A very fast reaction kinetics was also observed in the initial stages in the presence of oxygen. Our results show that the gasification reactions under H2O and CO2 started at lower temperatures than O2 thereby creating a partially oxidized structure. Due to the presence of a large number of activation sites, very high rates of gasification were achieved with oxygen. These findings could help identify optimal oxidant compositions towards maximizing carbon gasification and minimizing CO2 emissions.

Molecular analysis of hydrogen-bond structures in polymer electrolyte membrane in polymer electrolyte fuel cells below freezing temperatures

Fuel cells that generate electricity by using hydrogen and oxygen have been attracted to achieve carbon neutralization in the world because fuel cells don’t produce the greenhouse gases like carbon dioxide. The Polymer Electrolyte Fuel Cell (PEFC) has been used for many applications, such as cars and home stationary power supplies. Protons are transported in the polymer electrolyte membrane (PEM) and the transportation efficiency affects power generation efficiency. The internal state of the PEM has a great impact on the proton transport; however, the state is still not unclear below freezing temperatures because contained water might be frozen. Below freezing temperatures, proton-hopping, which is the mediation of proton conduction by water molecules, is important for proton transport. Proton-hopping is the mediation of proton conduction by water molecules. Therefore, we aim to understand the internal state of PEM below freezing temperatures. We performed reactive force-field molecular dynamics (ReaxFF MD) simulations to incorporate the proton-hopping. Our simulations consisted of two parts. First, bulk water was simulated at different temperatures to find the melting point of water by ReaxFF MD simulation. Radial distribution function (RDF) was used to find out whether the water is frozen or not. RDF of bulk water have a peak. When water is frozen, RDF of water have two distinct peaks and water molecules do not have translational or rotational motion. The size of the simulation box was 30.0 Å×30.0 Å×30.0 Å and the number of water molecules is 865. The temperature of the system was set at 1 K, 100 K, 200 K and 300 K. The time step was set at 0.50 fs and the sampling interval was 0.50 fs and the sampling interval was 100 steps. Second, the water condition in Nafion membrane was simulated. We used 4 Nafion polymer chains, 40 hydronium ions and 120 water molecules (water content: λ=4). The time step was set at 0.25 fs and the sampling interval was 100 steps. System of hydronium ions, water and Nafion changes faster than bulk water system and requires shorter timesteps. The temperature of the system was set at 300 K and 200 K. The annealing process was applied in order to equilibrate the system. After annealing, the production run was 250 ps using NVT ensemble. We performed the same analysis of the first simulation, and we analyzed the condition of water in the Nafion membrane. The RDF from the first simulation is shown in left figure. This figure shows that the two peaks were clearly visible at 1 K, 100 K and 200 K, although the second peak was barely visible at 300 K. Based on these results, the melting point of water in the ReaxFF potential was expected between 200 K and 300 K. The results of the second simulation are shown on the right figure. The results show that the heights of the two peaks are higher at 200 K than at 300 K. From these results, we could not find out whether the water in Nafion membrane was frozen or not only by these analysis methods. In this study, the melting point of water in the ReaxFF potential was found, and the state of water in polymer below freezing temperatures was analyzed. It was discussed in a paper by Qiangu Yan et al. According to them, water inside the Nafion membrane can be divided into three states. These are free water, bound water, and antifreeze water. While free water freezes, bound water has a lower melting point than normal water, and antifreeze water is not freeze. The specific heat of water in polymers should be analyzed in the future. Proton diffusion considering proton-hopping should also be analyzed in the future. Figure 1

Reactive molecular dynamics simulations of plastics pyrolysis with additives: Comparison of ReaxFF branches and experimental results

In this work, the efficiency of the combustion and aqueous branches of the reactive potential ReaxFF are compared through pyrolysis investigations from reactive molecular simulations of three polymers: polyethylene (PE), polystyrene (PS) and polypropylene (PP) with additive compounds. Thermogravimetric analyses were implemented, computing the weight loss percentage along with analyzes of the pyrolysis products. Three force fields were considered: the CHON-2019 and the CHOSFClN potentials coming from the combustion branch and the CHONCl-2022_weak potential coming from the aqueous branch. The results show that both CHON-2019 and CHONCl-2022_weak force fields give similar results in adequacy with the experimental results and are better than the CHOSFClN force field from the combustion branch. They can correctly predict the kinetic order of polymer degradation and identify the species formed and their quantities. As a result, the latest versions of the ReaxFF potential, belonging to the aqueous branch, are well suited to study the pyrolysis processes and may give equal or better results than force fields from the combustion branch. Furthermore, force fields from the aqueous branch seem to be more reliable for the investigation of the effect of the additives’ role on plastic pyrolysis. Finally, an important aspect presented in this article is the use of statistical simulations for the different investigations. An uncertainty estimation of 30 K was computed from ten thermogravimetric analyses simulated with ReaxFF on the PP pyrolysis.

Molecular dynamics simulation of deposition of amorphous carbon films on sapphire surfaces

The growth of amorphous carbon films on a sapphire surface was investigated using classical molecular dynamics simulation. The kinetic energy of carbon particles was set as 10 eV and ReaxFF potential was used to express the interaction between different kinds of particles. The results of the temperature distribution in both deposition time and deposition space are reported. Simulation results reveal that the grown amorphous carbon film consists of four regions, namely interlayer, low density, stable growth, and surface regions. In the interlayer region, the interlayer between substrate and pure carbon film is formed. In the low density region, a pure carbon film is grown while the film density decreases initially and then increases. In the stable growth region, the film density remains almost constant. The film density decreases rapidly in the surface region. The radial distribution function (RDF) analysis suggests that a structure similar to that of diamond exists in the stable growth region of the film. The lower film density in the low density and surface regions was interpreted to indicate the existence of abundant sp1 chain structures, which is supported by the depth profile of the sp fractions. The present results are in good agreement with previous experimental and simulation results and demonstrate the suitability of the ReaxFF potential in the simulation of amorphous carbon growth on sapphire substrate. Our study provides a good starting point for the simulation study of amorphous carbon films on sapphire substrates.

Ultralow-energy amorphization of contaminated silicon samples investigated by molecular dynamics.

Ion beam processes related to focused ion beam milling, surface patterning, and secondary ion mass spectrometry require precision and control. Quality and cleanliness of the sample are also crucial factors. Furthermore, several domains of nanotechnology and industry use nanoscaled samples that need to be controlled to an extreme level of precision. To reduce the irradiation-induced damage and to limit the interactions of the ions with the sample, low-energy ion beams are used because of their low implantation depths. Yet, low-energy ion beams come with a variety of challenges. When such low energies are used, the residual gas molecules in the instrument chamber can adsorb on the sample surface and impact the ion beam processes. In this paper we pursue an investigation on the effects of the most common contaminant, water, sputtered by ultralow-energy ion beams, ranging from 50 to 500 eV and covering the full range of incidence angles, using molecular dynamics simulations with the ReaxFF potential. We show that the expected sputtering yield trends are maintained down to the lowest sputtering yields. A region of interest with low damage is obtained for incidence angles around 60° to 75°. We also demonstrate that higher energies induce a larger removal of the water contaminant and, at the same time, induce an increased amorphization, which leads to a trade-off between sample cleanliness and damage.

Effect of particle size on ignition and oxidation of single aluminum: molecular dynamics study

Alumina nanoparticle is one of the attractive nanoparticles synthesized by the plasma method. The oxidation step in this method is challenging to explain experimentally. This work was to perform a molecular dynamics simulation to determine the oxidation mechanism of aluminum nanoparticles with different sizes and oxidation levels in the oxide layer. This work was to perform a molecular dynamics simulation to determine the oxidation mechanism of aluminum nanoparticles with different sizes and oxidation levels in the oxide layer. The simulation method employed the ReaxFF potential. The material used is aluminum nanoparticles in three different sizes (8, 12, and 16 nm) with an oxide layer thickness of 0.5 nm. Aluminum nanoparticles were given a relaxation treatment of 300 K for 1 ps and then heated to a temperature of 3250 K with a heating rate of 5×1013 K/s and cooled to 300 K. The ensemble used is a canonical ensemble with the Nose/Hoover thermostat method. The result shows that the higher the temperature applied to the system, the more oxygen molecules adsorption occurs on the surface of the oxide layer and the diffusion of oxygen to the particle core. The higher temperature applied also causes gaps, or void spaces, between the core and the shell. The reaction barrier for diffusion of oxygen also decreased significantly due to void space, and the surface of the aluminum core dissociates to the surface (alumina shell). Particles with a smaller size have a shorter ignition delay time. In addition, the smaller the particle size, the more oxygen molecules' reacted with aluminum particles in the particle core

Molecular Dynamics Modeling of SiO2 Melts and Glass Formation Processes

Molecular dynamics (MD) with ReaxFF potentials is used to study the melting process of quartz and cristobalite together with the amorphous structures obtained at different stages of melting by cooling the melt. The long-term preservation of an excess of eight-membered rings inherited from the crystalline phase is found in the quartz melts, while in the cristobalite melts, the similar preservation of six-membered rings is not observed. Thus, it can be stated that the quartz melts and glasses obtained from them have structural memory, in contrast to cristobalite melts. An increase in the number of four-membered rings with increasing temperature is revealed. A number of other features of the obtained amorphous structures, which we consider as models for glasses, are discussed.

The Dynamic Nature of Graphene Active Sites in the H2O Gasification process: A ReaxFF and DFT Study

A steam-rich environment is a more promising application scenario for future coal-fired processes, while active sites are the key factor that determines the reactivity of carbonaceous fuels. The steam gasification process of carbon surfaces with different numbers of active sites (0, 12, 24, 36) was simulated using reactive molecular dynamics in the present study. The temperature for the decomposition of H2O and the gasification of carbon is determined using temperature-increasing simulation. The decomposition of H2O was influenced by two driving forces, thermodynamics and active sites on the carbon surface, which dominated the different reaction stages, leading to the observed segmentation phenomenon of the H2 production rate. The existence and number of initial active sites have a positive correlation with both two stages of the reaction, greatly reducing the activation energy. Residual OH groups play an important role in the gasification of carbon surfaces. The supply of OH groups through the cleavage of OH bonds in H2O is the rate-limiting step in the carbon gasification reaction. The adsorption preference at carbon defect sites was calculated using density functional theory. Two stable configurations (ether & semiquinone groups) can be formed with O atoms adsorbed on the carbon surface according to the number of active sites. This study will provide further insights into the tuning of active sites for advanced carbonaceous fuels or materials. The large-scale atomic/molecule massively parallel simulator (LAMMPS) code combined with the reaction force-field method was used to carry out the ReaxFF molecular dynamics simulation, where the ReaxFF potentials were taken from Castro-Marcano, Weismiller and William. The initial configuration was built using Packmol, and the visualization of the calculation results was realized through Visual Molecular Dynamics (VMD). The timestep was set to 0.1 fs to detect the oxidation process with high precision. PWscf code in QUANTUM ESPRESSO (QE) package, was used to evaluate the relative stability of different possible intermediate configurations and the thermodynamic stability of gasification reactions. The projector augmented wave (PAW) and the generalized gradient approximation of Perdew-Burke-Ernzerhof (PBE-GGA) were adopted. Kinetic energy cutoffs of 50 Ry and 600 Ry, and a uniform mesh of 4 × 4 × 1 k-points were used.

Density Functional Tight-Binding Model for Lithium-Silicon Alloys.

The predictive power of molecular dynamic simulations is mainly restricted by the time scale and model accuracy. Many systems of current relevance are of such complexity that they require addressing both issues simultaneously. This is the case of silicon electrodes in Li-ion batteries, where different LixSi alloys are formed during charge/discharge cycles. While first-principles treatments for this system are seriously limited by the computational cost of exploring its large conformational space, classical force fields are not transferable enough to represent it accurately. Density Functional Tight Binding (DFTB) is an intermediate complexity approach capable of capturing the electronic nature of different environments with a relatively low computational cost. In this work, we present a new set of DFTB parameters suited to model amorphous LixSi alloys. LixSi is the usual finding upon cycling the Si electrodes in the presence of Li ions. The model parameters are constructed with a particular emphasis on their transferability for the entire LixSi composition range. This is achieved by introducing a new optimization procedure that weights stoichiometries differently to improve the prediction of their formation energies. The resulting model is shown to be robust for predicting crystal and amorphous structures for the different compositions, giving excellent agreement with DFT calculations and outperforming state-of-the-art ReaxFF potentials.