kMC Simulations Research Articles

Multiscale modeling of damaged surface topology in a hypersonic boundary.

In this work, we used molecular dynamics (MD) to perform trajectory simulations of ice-like argon and amorphous silica aggregates on atomically smooth highly ordered pyrolytic graphite (HOPG) and a comparatively rougher quartz surface. It was found that at all incidence velocities, the quartz surface was stickier than the HOPG surface. The sticking probabilities and elastic moduli obtained from MD were then used to model surface evolution at a micron length scale using kinetic Monte Carlo (kMC) simulations. Rules were derived to control the number of sites available for the process execution in kMC to accurately model erosion of HOPG by atomic oxygen (AO) attack and ice-nucleation on surfaces. It was observed that the effect of defects was to increase the material erosion rate, while that of aggregate nucleation was to lower it. Similarly, simulations were performed to study the effects of AO attack and N2 adsorption-desorption on surface evolution and it was found that N2 adsorption-desorption limits the surface available for erosion by AO attack.

Efficient Implementation of Cluster Expansion Models in Surface Kinetic Monte Carlo Simulations with Lateral Interactions: Subtraction Schemes, Supersites, and the Supercluster Contraction.

While lateral interaction models for reactions at surfaces have steadily gained popularity and grown in terms of complexity, their use in chemical kinetics has been impeded by the low performance of current kinetic Monte Carlo (KMC) algorithms. The origins of the additional computational cost in KMC simulations with lateral interactions are traced back to the more elaborate cluster expansion Hamiltonian, the more extensive rate updating, and to the impracticality of rate-catalog-based algorithms for interacting adsorbate systems. Favoring instead site-based algorithms, we propose three ways to reduce the cost of KMC simulations: (1) representing the lattice energy by a smaller Supercluster Hamiltonian without loss of accuracy, (2) employing the subtraction schemes for updating key quantities in the simulation that undergo only small, local changes during a reaction event, and (3) applying efficient search algorithms from a set of established methods (supersite approach). The cost of the resulting algorithm is fixed with respect to the number of lattice sites for practical lattice sizes and scales with the square of the range of lateral interactions. The overall added cost of including a complex lateral interaction model amounts to less than a factor 3. Practical issues in implementation due to finite numerical accuracy are discussed in detail, and further suggestions for treating long-range lateral interactions are made. We conclude that, while KMC simulations with complex lateral interaction models are challenging, these challenges can be overcome by modifying the established variable step-size method by employing the supercluster, subtraction, and supersite algorithms (SSS-VSSM). © 2019 Wiley Periodicals, Inc.

Parallel BFS implementing optimized decomposition of space and kMC simulations for diffusion of vacancies for quantum storage

We present a novel breadth-first search (BFS) algorithm based on the notion of temporal evolvability that is adaptable to various multicore architectures for simulating diffusion of vacancies in hexagonal silicon carbide (4H-SiC) for information storage. The algorithm is formulated in the semi-ring algebraic framework of BFS and incorporates a real-space grid decomposition to optimize the number of nodes that are evaluated in each frontier of the evolution.Scaling characteristics are first evaluated from performance runs for two formulations: (i) recursive depth-first search (DFS) and (ii) semi-ring implementations of BFS. The results for a real-space grid implementation of BFS are then presented. We demonstrate that each new iteration reduces the communication overhead, enhancing performance. A comparison of the real-space grid based BFS with a kinetic Monte Carlo (kMC) algorithm is also presented for the case of diffusion without the influence of the Coulomb interaction. The efficient parallel implementations of this latter approach enables simulations of larger systems than those studied before.

Charge Pair Separation Dynamics in Organic Bulk‐Heterojunction Solar Cells

Charge Pair Separation Dynamics in Organic Bulk‐Heterojunction Solar Cells

Characterising submonolayer deposition via visibility graphs

We use visibility graphs as a tool to analyse the results of kinetic Monte Carlo (kMC) simulations of submonolayer deposition in a one-dimensional point island model. We introduce an efficient algorithm for the computation of the visibility graph resulting from a kMC simulation and show that from the properties of the visibility graph one can determine the critical island size, thus demonstrating that the visibility graph approach, which implicitly combines size and spatial data, can provide insights into island nucleation and growth processes.

Kinetic Monte Carlo simulations of coverage effect on Ag and Au monolayers growth on Cu (1 1 0)

In this work, different coverages (1% ≤ coverage ≤ 20%) are studied during the growth of Ag and Au monolayers on the Cu (1 1 0) surface using kinetic Monte Carlo simulations (kMC). The results have shown that as coverage increases, the likelihood of an isolated adatom connecting to another adatom or island also increases. Thus, the phenomenon by which an adatom disintegrates from one island to another becomes significant. Molecular statics method is used to calculate the activation energy barriers of several diffusion processes in the framework of the Embedded-Atom method (EAM). These diffusion processes are implemented in the kinetic Monte Carlo (kMC) simulations under the Anisotropic Bond-Breaking Model (ABBM) model. The results show that as the coverage increases, 1D clusters are observed at 200 K. However, the transition from compact 1D (at 200 K) to 2D islands has been observed for relatively high temperatures (300 K and 400 K).

Highly Efficient Ni-Doped Iron Catalyst for Ammonia Synthesis from Quantum-Mechanics-Based Hierarchical High-Throughput Catalyst Screening

To discover more efficient industrial catalysts for ammonia synthesis via the Haber–Bosch (HB) process, we employed quantum mechanics (QM)-based hierarchical high-throughput catalyst screening (HHTCS) to test a wide group of elements (34) as candidates to dope the Fe(111) catalyst subsurface. The QM free-energy reaction network of HB over Fe(111) yields ten barriers as potentially rate-determining, of which we select four as prototypical, arrange them hierarchically, and define a corresponding set of screening criteria, which we then use to screen candidate catalysts. This leads to two promising candidates (Co and Ni), from which we selected the most promising (Ni) for a complete QM and kinetic study. The kinetic Monte Carlo (kMC) simulations predict a 16-fold increase in HB turn-over frequency (TOF) for the Ni-doped catalyst compared to the pure Fe(111) surface under realistic conditions. The 16-fold increase in HB TOF is a significant improvement and may trigger future experimental studies to validate our prediction. This TOF improvement could lead to similar reaction rates as with pure Fe but at a reaction temperature decreased by 100° from 773 to 673 K and a total reactant pressure decreased by 6 times from 201 to 34 atm. We interpret the reasons underlying this improvement using valence bond and kinetic analyses. We suggest this Ni-doped Fe(111) catalyst as a candidate to reduce the world energy consumption for the HB process while satisfying future needs for energy and environment.

Hierarchical analysis of alloying element effects on gas nitriding rate of Fe alloys: A DFT, microkinetic and kMC study

Nitriding is the most widely employed thermochemical surface treatment to enhance the mechanical properties of steel. Specifically, gas nitriding, which is a low-temperature process for efficiently producing high-performance steels, has a disadvantage in that it consumes a large amount of time. To enhance the nitriding rate, we studied the surface alloying of iron (Fe) and its effect on ammonia (NH3) nitriding of Fe using a hierarchical protocol with density functional theory (DFT)-based microkinetics and real-time simulations. First, we considered the NH3 decomposition and nitrogen (N) diffusion mechanism on clean and alloyed (Fe-X) Fe (100) surfaces using DFT. In this study, the alloying elements including transition metals and period III to VI elements in the periodic table were considered for DFT-based computational screening. For the candidate Fe-X systems selected to improve the nitriding rate in the previous step, we calculated all the energy barriers for every elementary reaction step by varying the alloying elements and performed microkinetic analysis using those kinetic energy barriers to determine their influence on the nitriding rate. After adding consideration of thermodynamic factors, selected candidate alloys were subjected to detailed DFT calculations of the nitriding mechanism with N coverage, and based on these results, a kinetic Monte Carlo (kMC) simulation was performed to reconfirm the results under the actual nitriding process conditions. Through a hierarchical protocol, we performed a theoretical analysis and simulation of the effects of alloying elements on the nitriding rate that were not explained experimentally and suggested the best alloying element with the improved nitriding rate.

Dynamics of 3D-island growth on weakly-interacting substrates

The growth dynamics of faceted three-dimensional (3D) Ag islands on weakly-interacting substrates are investigated—using kinetic Monte Carlo (kMC) simulations and analytical modelling—with the objective of determining the critical top-layer radius Rc required to nucleate a new island layer as a function of temperature T, at a constant deposition rate. kMC shows that Rc decreases from 17.3 to 6.0 Å as T is increased at 25 K intervals, from 300 to 500 K. That is, a higher T promotes top-layer nucleation resulting in an increase in island height-to-radius aspect ratios. This explains experimental observations for film growth on weakly-interacting substrates, which are not consistent with classical homoepitaxial growth theory. In the latter case, higher temperatures yield lower top-layer nucleation rates and lead to a decrease in island aspect ratios. The kMC simulation results are corroborated by an analytical mean field model, in which Rc is estimated by calculating the steady-state adatom density on the island side facets and top layer as a function of T. The overall findings of this study constitute a first step toward developing rigorous theoretical models, which can be used to guide synthesis of metal nanostructures, and layers with controlled shape and morphology, on technologically important substrates, including two-dimensional crystals, for nanoelectronic and catalytic applications.

Multiscale modelling from quantum level to reactor scale: An example of ethylene epoxidation on silver catalysts

Abstract Ethylene epoxidation is one of the most important selective chemical oxidations in industry. For a controlled transformation of ethylene (ethene) into epoxide, silver is the only commercially suitable catalyst. Although it is usually doped, even with pristine silver activity, selectivity, and stability vary strongly with facets. In this work, we use this reaction on Ag( 111 ) and Ag( 100 ) as a classical formation model to demonstrate the capabilities of physical multiscale modelling, to show why Ag( 100 ) nanocubes offer superior catalysis, and to optimise reactivity. First, we describe the elementary reactions on pristine surfaces with the quantum chemistry calculations, using density functional theory (DFT). The free energies of all intermediates, kinetic rates from the transition state theory and adsorption/desorption equilibria are calculated from first principles. These results are applied to kinetic Monte Carlo (kMC) simulations, where the spatio-temporal evolution of the system on a meso-scale can be followed. The differences in activity, concentration, selectivity, and apparent activation energy are observed, investigated, and analysed. Lastly, mean-field concepts – micro-kinetics and computational fluid dynamics (CFD) – are used to simulate how the synthesis proceeds in a reactor. Mechanism, catalytic coverage and the effects of pressure, temperature, and particle composition, size and shape on the performance are evaluated. We show that multiscale modelling is a powerful instrumental approach for real unit engineering, while the level of detail required is dictated by the purpose of a representation and available resources.

Hierarchical Analysis of Alloying Element Effects on Gas Nitriding Rate of Fe Alloys: A DFT, Microkinetic and kMC Study

Nitriding is the most widely employed thermochemical surface treatment to enhance the mechanical properties of steel. Specifically, gas nitriding, which is a low-temperature process for efficiently producing high-performance steels, has a disadvantage in that it consumes a large amount of time. To enhance the nitriding rate, we studied the surface alloying of iron (Fe) and its effect on ammonia nitriding of Fe using a hierarchical protocol with density functional theory (DFT)-based microkinetics and real-time simulations. First, we considered the NH3 decomposition and nitrogen (N) diffusion mechanism on clean and alloyed (Fe-X) Fe (100) surfaces using DFT. In this study, the alloying elements including transition metals and period III to VI elements in the periodic table were considered for DFT-based computational screening. For the candidate Fe-X systems selected to improve the nitriding rate in the previous step, we calculated all the energy barriers for every elementary reaction step by varying the alloying elements and performed microkinetic analysis using those kinetic energy barriers to determine their influence on the nitriding rate. After adding consideration of thermodynamic factors, selected candidate alloys were subjected to detailed DFT calculations of the nitriding mechanism with N coverage, and based on these results, a kinetic Monte Carlo (kMC) simulation was performed to reconfirm the results under the actual nitriding process conditions. Through a hierarchical protocol, we performed a theoretical analysis and simulation of the effects of alloying elements on the nitriding rate that were not explained experimentally and suggested the best alloying element with the improved nitriding rate.

Microscopic modeling and optimal operation of thermal atomic layer deposition

This work develops a comprehensive framework for first-principles-based microscopic modeling, data-driven input/output modeling and optimal operation of thermal atomic layer deposition (ALD) of SiO2 thin-films using bis(tertiary-butylamino)silane (BTBAS) and ozone as precursors. Specifically, we initially utilize Density Functional Theory (DFT)-based calculations for the computation of the key thermodynamic and kinetic parameters, which are then used in the microscopic modeling of the ALD process. Subsequently, a detailed microscopic model is constructed, accounting for the microscopic lattice structure and atomic interactions, as well as multiple microscopic film growth processes including physisorption, abstraction and competing chemical reaction pathways. Kinetic Monte-Carlo (kMC) algorithms are utilized to obtain computationally efficient microscopic model solutions while preserving model fidelity. The obtained kMC simulation results are used to train Artificial Neural Network (ANN)-based data-driven models that capture the relationship between operating process parameters and time to ALD cycle completion. Specifically, a dense two-hidden-layer feed-forward ANN is constructed to find a feasible range of ALD operating conditions accounting for industrial considerations, and a Bayesian Regularized ANN is constructed to implement the cycle-to-cycle optimization of ALD cycle time. Extensive simulation results demonstrate the effectiveness of the proposed approaches. The kMC models successfully achieves a growth per cycle (GPC) of 1.8Å per cycle, which is in the range of reported experimental values. The ANN models accurately predict deposition time to steady-state from the given operating condition input, and the cycle time optimization using BRANN model reduces the conventional BTBAS cycle time by 60%.

Effect of Co doping on mechanism and kinetics of ammonia synthesis on Fe(1 1 1) surface

With the aim of improving the efficiency of the Haber-Bosch industrial process for the synthesis of ammonia, here we explore doping the traditional Fe-based Haber-Bosch catalyst with an impurity element. Starting from a previous experimentally-validated theoretical investigation of the reaction mechanism for Haber-Bosch synthesis of ammonia on the Fe bcc(1 1 1) surface, we focus on changes in mechanism and kinetics brought about by substitutional doping of 25% top layer iron with cobalt. The choice of Co is justified by the analysis of the wave functions of the critical reaction steps on the Fe(1 1 1) surface which showed that large changes in the net spin (magnetization) of the Fe atoms are thereby involved, and suggested that dopants with modified spins might accelerate rates. Quantum Mechanics values of free energies and reaction barriers are calculated for the Co-doped system for a set of 20 important surface configurations of adsorbates, and used as input to kinetic Monte Carlo (kMC) simulations to obtain final ammonia production. We find that at T = 673 K, P(H2) = 15 atm, P(N2) = 5 atm, and P(NH3) = 1 atm, target conditions to drastically reduce the extreme energy cost of industrial ammonia synthesis process, top-layer Co doping leads to an acceleration by a factor of 2.3 in reaction rates of ammonia synthesis, and therefore an expected corresponding decrease in production costs.

Relaxation in one-dimensional chains of interacting magnetic nanoparticles: Analytical formula and kinetic Monte Carlo simulations

We study the relaxation of long chains of magnetic nanoparticles (MNPs). In spite of the simplicity of this system, there is no theoretical framework for this basic assembly. Using the two-level approximation for energy, we perform first-principles calculations and kinetic Monte Carlo (kMC) simulations to obtain the effective relaxation time ${\ensuremath{\tau}}_{N}^{}$ of the chain by incorporating the effects of dipole-dipole interactions and anisotropy axes orientation of the MNPs. For analytical tractability, we consider the case when all easy-axis and initial magnetic moments make an angle $\ensuremath{\alpha}$ with the chain axis. In the absence of dipolar interactions, the relaxation is governed, as expected, by the usual N\'eel relaxation time ${\ensuremath{\tau}}_{N}^{0}$. In the presence of interactions, the magnetic relaxation curve is always perfectly fitted by an exponentially decaying function. The dipolar field induces antiferromagnetic or ferromagnetic interactions between the moments: depending on $\ensuremath{\alpha}$ values, this induces a fastening of relaxation time $\left({\ensuremath{\tau}}_{N}^{}l{\ensuremath{\tau}}_{N}^{0}\right)$ or a slowing down $\left({\ensuremath{\tau}}_{N}^{}g{\ensuremath{\tau}}_{N}^{0}\right)$. The analytical determination of ${\ensuremath{\tau}}_{N}^{}$ is nontrivial, but we have obtained an approximate form that is confirmed by kMC simulations. Finally, it is shown that the equilibrium state is comprised of short-lived ferromagnetic and antiferromagnetic domains, the size of which increases with the dipolar strength. We believe that the above conclusions can be drawn for chains with more complicated structures exhibiting bends, curls, and intersections in higher dimensions. Our study is relevant in the context of applications such as magnetic recording, digital data processing, and magnetic hyperthermia, in which long chains of MNPs are ubiquitous.

Morphology evolution of magnesium facets: DFT and KMC simulations

One of the crucial steps for the development of batteries is understanding the interface stability and morphological changes occurring during continuous stripping and deposition. In order to investigate the dependence of morphology evolution on surface orientation, we examine the energetics and growth mechanism on magnesium (0001), (101[combining macron]0), (101[combining macron]1), (112[combining macron]0) and (112[combining macron]1) surface orientations using density functional theory and kinetic Monte Carlo simulations. Workfunctions, surface, adsorption and interaction energies, diffusion barriers and k-rates for diffusion via hopping and exchange mechanisms are studied. The results provide a comprehensive relationship between these properties and morphology evolution. The latter shows strong dependence on the surface orientation, demonstrating the need for all commonly present facets to be studied, instead of focusing only on the most stable surface orientation.

Phase Transformation Assisted Twinning in Face-Centered-Cubic FeCrNiCoAl <sub>0.36</sub>

The FeNiCoCr-based high entropy alloys (HEAs) exhibit excellent mechanical properties, such as twin-induced plasticity (TWIP) and phase transformation plasticity (TRIP) that can reach a remarkable combination of strength and ductility. In this work, the face-centered-cubic (FCC) single-crystal FeNiCoCrAl0.36 HEAs were studied, using the density functional theory (DFT) combined with the phonon calculation to estimate the stacking fault energies, temperature-dependent phase stabilities of different structures. And the kinetic Monte Carlo (kMC) is used to predict the substructures evolution based on the transition state energies obtained from DFT calculations. We proposed two different formation paths of nano-twins in this Al-composited HEA and found that short-range hexagonal-close-packed (HCP)-stacking could occur in this HEA. DFT calculations suggest that this HEA has negative stacking fault energy, HCP formation energy, and twin-formation energy at 0 K. Phonon calculations represent that at the finite temperature, the competing FCC/HCP phase stability and propensity for twinning makes it possible to form HCP-like twin boundaries. The kMC simulations suggest that under deformation, HCP substructures could form followed by twins which differs to the study of others. With the great agreement of results from kMC simulations and experiments, this twin-formation path offers a new concept of designing TWIP HEAs containing tunable twin structures with HCP and TWIN lamellae structures, which results in better mechanical properties of HEAs.

Complex oscillatory decrease with size in diffusivity of {100}-epitaxially supported 3D fcc metal nanoclusters.

Diffusion and coalescence of supported 3D metal nanoclusters (NCs) leads to Smoluchowski Ripening (SR), a key pathway for catalyst degradation. Variation of the NC diffusion coefficient, DN, with size N (in atoms) controls SR kinetics. Traditionally, a form DN∼N-β was assumed consistent with mean-field analysis. However, KMC simulation of a stochastic model for diffusion of {100}-epitaxially supported fcc NCs mediated by surface diffusion reveals instead a complex oscillatory decrease of DN with N. Barriers for surface diffusion of metal atoms across and between facets, along step edges, etc., in this model are selected to accurately capture behavior for fcc metals. (This contrasts standard bond-breaking prescriptions which fail dramatically.) For strong adhesion, equilibrated NCs are truncated pyramids (TP). Local minima of DN sometimes but not always correspond to sizes, NTP, where these have a closed-shell structure. Local maxima generally correspond to N≈NTP + 3 for N = O(102). For weak adhesion, equilibrated NCs are truncated octahedra (TO), and local minima of DN occur for sizes close or equal to those of just a subset of closed-shell structures. Analytic characterization of energetics along the NC diffusion pathway (which involves dissolving and reforming outer layers of facets) provides fundamental insight into the behavior of DN, including the strong variation with N of the effective NC diffusion barrier.

Direct coupling of continuum and kinetic Monte Carlo models for multiscale simulation of electrochemical systems

Electrochemical systems include atomistic processes at electrochemical interfaces and macroscopic transport processes, which can be modeled using kinetic Monte Carlo (kMC) simulation and continuum equations, respectively. Multiparadigm algorithms are applied to directly couple such models to study multiscale interactions. This article compares different algorithms for an example problem. Results quantify the effect of computational cost and numerical accuracy by the choice of algorithm and its configuration. The stochastic fluctuations of kMC simulations as well as sequential data exchange between the models generate errors in coupled simulations. Measures to reduce stochastic fluctuations or revise exchanged data can be either highly successful or futile, depending on the dominant cause of the error. Hence, we strongly advise to identify the different causes of errors and their mechanics when selecting a coupling algorithm or optimizing its configuration. This article provides various algorithms and suggestions for their configuration to enable efficient and robust multiscale simulations.

Machine Learning–Based Charge Transport Computation for Pentacene

AbstractInsight into the relation between morphology and transport properties of organic semiconductors can be gained using multiscale simulations. Since computing electronic properties, such as the intermolecular transfer integral, using quantum chemical (QC) methods requires a high computational cost, existing models assume several approximations. A machine learning (ML)–based multiscale approach is presented that allows to simulate charge transport in organic semiconductors considering the static disorder within disordered crystals. By mapping fingerprints of dimers to their respective transfer integral, a kernel ridge regression ML algorithm for the prediction of charge transfer integrals is trained and evaluated. Since QC calculations of the electronic structure must be performed only once, the use of ML reduces the computation time radically, while maintaining the prediction error small. Transfer integrals predicted by ML are utilized for the computation of charge carrier mobilities using off‐lattice kinetic Monte Carlo (kMC) simulations. Benefiting from the rapid performance of ML, microscopic processes can be described accurately without the need for phenomenological approximations. The multiscale system is tested with the well‐known molecular semiconductor pentacene. The presented methodology allows reproducing the experimentally observed anisotropy of the mobility and enables a fast estimation of the impact of disorder.

Adlayer structure and lattice size effects on catalytic rates predicted from KMC simulations: NO oxidation on Pt(111).

Repulsive and/or attractive interactions between surface adsorbates have an important effect on the structure of the adsorbate layer and consequently on the rate of heterogeneous catalytic reactions. Thus, developing reaction models that take into account adsorbate-adsorbate interactions is crucial for making accurate predictions of the catalytic rate and surface coverage during reaction. In the present work, we employ kinetic Monte Carlo simulation to model the catalytic NO oxidation on Pt (111), adopting a cluster expansion (CE) Hamiltonian approach for treating the aforementioned interactions. We investigate CEs of increasing complexity, ranging from pairwise 1st nearest neighbor to long-range and many-body terms. We show that energetic models incorporating solely short-range interactions result in ordered adlayer structures, which are disrupted by anti-phase boundaries and defective regions when the size of the periodic lattice is non-commensurate to the structure of the stable adlayer. We find that O2 dissociates on sites located in these defective regions, which are predominantly responsible for the activity, and the predicted catalytic rate is strongly depended on the lattice size. Such effects are absent when employing non-periodic lattices, whereon the catalytic activity appears more intense on edges/corner sites. Finally, inclusion of long-range interactions in the model Hamiltonian induces relative disorder in the adsorbate layer, which is ascribed to the "softening" of the repulsive interactions between adspecies. Under these circumstances, the distribution of activation energies for O2 dissociation is broader as compared to short-range interaction models and on this basis we explain the disparate catalytic rate predictions when using different CEs.