Cluster Expansion Hamiltonians Research Articles

KMCpy: A python package to simulate transport properties in solids with kinetic Monte Carlo

Understanding ion transport in functional materials is crucial to unravel complex chemical reactions, improving the rate performance of materials for energy storage and conversion, and optimizing catalysts. To model ion transport, atomistic simulations, including molecular dynamics (MD) and kinetic Monte Carlo (kMC) have been developed and applied. However, kMC simulations are utilized to a lower extent than MDs due to a lack of systematic workflows to construct models for predicting transition rates. Here, we present kMCpy, a lightweight, customizable, and modular python package to compute the ionic transport properties in crystalline materials using kMC. kMCpy is remarkably versatile and user-friendly, making it a powerful code for studying materials′ kinetics in crystalline systems. kMCpy can be combined with (local) cluster expansion Hamiltonians derived from first-principles calculations. kMCpy is versatile with respect to any type of crystalline material, bearing any dimensionality, i.e., 1D, 2D, and 3D. kMCpy provides (i) a comprehensive workflow to enumerate all possible migration events in crystalline systems, (ii) to derive transition rates efficiently and at the accuracy of first-principles calculations, and (iii) a robust kMC solver to study kinetic phenomena in materials. The workflow implemented in kMCpy provides a systematic way to compute highly accurate kinetic properties. Hence, kMCpy can be used in high-throughput simulations for the discovery and optimization of novel functional materials.

Revealing the energetic and kinetic stabilisation of Guinier–Preston zone morphology during growth in binary Al–Cu alloys

In this paper we are going to present a computational study on the growth behaviour of Guinier–Preston zones (GPZ) in binary Al–Cu alloys. For this purpose we performed atomistic Monte-Carlo (MC) simulations and calculations of static configurational energies by applying a Cluster Expansion Hamiltonian for the binary Al–Cu system derived from density functional theory, which we already developed in a former publication. The MC simulations were designed to represent Cu diffusion and growth events of GPZ and were utilised to display atomistically resolved concentration fields and diffusion currents of mobile Cu in the vicinity of GPZ. For the static calculations we constructed sequences of configurations mimicking growth paths of Cu atoms and recorded their formation energy profiles. From both parts an energetically as well as kinetically stabilised mechanism for GPZ growth is proposed and the origin for the transition from GP(I) to GP(II) is elucidated.

A Caching Scheme To Accelerate Kinetic Monte Carlo Simulations of Catalytic Reactions

Kinetic Monte Carlo (KMC) simulations have been instrumental in advancing our fundamental understanding of heterogeneously catalyzed reactions, with particular emphasis on structure sensitivity, ensemble effects, and the interplay between adlayer structure and adsorbate-adsorbate lateral interactions in shaping the observed kinetics. Yet, the computational cost of KMC remains high, thereby motivating the development of acceleration schemes that would improve the simulation efficiency. We present an exact such scheme, which implements a caching algorithm along with shared-memory parallelization to improve the computational performance of simulations incorporating long-range adsorbate-adsorbate lateral interactions. This scheme is based on caching information about the energetic interaction patterns associated with the products of each possible lattice process (adsorption, desorption, reaction etc.). Thus, every time a reaction occurs ("ongoing reaction"), it enables fast updates of the rate constants of "affected reactions", i.e., possible reactions in the region of influence of the "ongoing reaction". Benchmarks on KMC simulations of NOx oxidation/reduction, yielded acceleration factors of up to 20, when comparing single-thread runs without caching to runs on 16 threads with caching, for simulations with a cluster expansion Hamiltonian that incorporates up to 8th-nearest-neighbor interactions.

Role of short- and long-range ordering on diffusion in Ni-Al alloys

The varying degrees of short- and long-range order exhibited by multicomponent solids complicates a first-principles calculation of nondilute diffusion coefficients. Temporal and spatial variations in the local degree of order affect the migration barriers of individual hops and can result in strong correlations between successive hops that ultimately affect macroscopic transport coefficients. Here we report on a first-principles study of diffusion in Ni-rich Ni-Al alloys. We used cluster expansion Hamiltonians to describe the energies of the end states and migration barriers of each hop. Kinetic Monte Carlo simulations were performed to calculate macroscopic transport coefficients. Variations in the degree of ordering are shown to play a significant role in affecting diffusion coefficients. While Al has a higher mobility than Ni in the disordered Ni-rich fcc solid solution, it becomes significantly less mobile in the ordered $L{1}_{2}\phantom{\rule{4pt}{0ex}}{\ensuremath{\gamma}}^{\ensuremath{'}}$ phase due to strong thermodynamic tendencies that keep Al trapped to disconnected sublattice sites.

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.

First principles phase diagram calculation for the 2D TMD system [formula omitted

First principles phase diagram calculations, that included van der Waals interactions, were performed for the bulk transition metal dichalcogenide system (1−X)·WS2−(X)·WTe2. To obtain a converged phase diagram, a series of cluster expansion calculations were performed with increasing numbers of structural energies, (Nstr) up to Nstr=435, used to fit the cluster expansion Hamiltonian. All calculated formation energies are positive and all ground-state analyses predict that formation energies for supercells with 16 or fewer anion sites are positive; but when 150⪅Nstr⪅376, false ordered ground-states are predicted. With Nstr≥399, only a miscibility gap is predicted, but one with dramatic asymmetry opposite to what one expects from size-effect considerations; i.e. the calculations predict more solubility on the small-ion S-rich side of the diagram and less on the large-ion Te-rich side. This occurs because S-rich low-energy metastable ordered configurations have lower energies than their Te-rich counterparts which suggests that elastic relaxation effects are not dominant for the shape of the miscibility gap.

Resolving phase stability in the Ti-O binary with first-principles statistical mechanics methods

The Ti-O system consists of a multitude of stable and metastable oxides that are used in wide ranging applications. In this work we investigate phase stability in the Ti-O binary from first principles. We perform a systematic search for ground state structures as a function of oxygen concentration by considering oxygen-vacancy and/or titanium-vacancy orderings over four parent crystal structures: (i) hcp Ti, (ii) $\ensuremath{\omega}$-Ti, (iii) rocksalt, and (iv) hcp oxygen containing interstitial titanium. We explore phase stability at finite temperature using cluster expansion Hamiltonians and Monte Carlo simulations. The calculations predict a high oxygen solubility in hcp Ti and the stability of suboxide phases that undergo order-disorder transitions upon heating. Vacancy ordered rocksalt phases are also predicted at low temperature that disorder to form an extended solid solution at high temperatures. Predicted stable and metastable phase diagrams are qualitatively consistent with experimental observations, however, important discrepancies are revealed between first-principles density functional theory predictions of phase stability and the current understanding of phase stability in this system.

A machine learning approach to graph-theoretical cluster expansions of the energy of adsorbate layers.

The accurate description of the energy of adsorbate layers is crucial for the understanding of chemistry at interfaces. For heterogeneous catalysis, not only the interaction of the adsorbate with the surface but also the adsorbate-adsorbate lateral interactions significantly affect the activation energies of reactions. Modeling the interactions of the adsorbates with the catalyst surface and with each other can be efficiently achieved in the cluster expansion Hamiltonian formalism, which has recently been implemented in a graph-theoretical kinetic Monte Carlo (kMC) scheme to describe multi-dentate species. Automating the development of the cluster expansion Hamiltonians for catalytic systems is challenging and requires the mapping of adsorbate configurations for extended adsorbates onto a graphical lattice. The current work adopts machine learning methods to reach this goal. Clusters are automatically detected based on formalized, but intuitive chemical concepts. The corresponding energy coefficients for the cluster expansion are calculated by an inversion scheme. The potential of this method is demonstrated for the example of ethylene adsorption on Pd(111), for which we propose several expansions, depending on the graphical lattice. It turns out that for this system, the best description is obtained as a combination of single molecule patterns and a few coupling terms accounting for lateral interactions.

Short-Range Order in High Entropy Alloys: Theoretical Formulation and Application to Mo-Nb-Ta-V-W System

In high-entropy alloys (HEAs), the local chemical fluctuations from disordered solute solution state into segregation, precipitation and ordering configurations are complex due to the large number of elements. In this work, the cluster expansion (CE) Hamiltonian for multi-component alloy systems is developed in order to investigate the dependence of chemical ordering of HEAs as a function of temperature dependence due to derivation of configuration entropy from the ideal solute solution. Analytic expressions for Warren–Cowley short-range order (SRO) parameters are derived for a five component alloy system. The theoretical formulation is used to investigate the evolution of the ten different SRO parameters in the MoNbTaVW and the sub-quaternary systems obtained by Monte-Carlo simulations within the combined CE and first-principles formalism. The strongest chemical SRO parameter is predicted for the first nearest-neighbor Mo-Ta pair that is in consistent agreement with high value of enthalpy of mixing in the B2 structure for this binary system. The prediction of B2 phase presence for Mo-Ta pairs in the considered bcc HEAs is reinforced by the positive contribution to the average SRO from the second nearest-neighbor shell. Interestingly, it is found that the average SRO parameter for the first and second nearest-neighbor shells of V-W pairs is also strongly negative in a comparison with the Mo-Ta pairs. This finding in the HEAs can be rationalized and discussed by the presence of the ordered-like B32 phase which has been predicted as the ground-state structure in binary bcc V-W system at the equimolar composition.

First-principles calculations of the phase equilibrium of BexZn1−xO alloys

Bandgap engineering of ZnO is crucial towards its practical applications. Due to their wide bandgap, BexZn1−xO alloys are promising materials for making optoelectronic devices that function in the solar-blind wavelength region. In this study, a theoretical investigation of the thermodynamics of these BexZn1−xO alloys is carried out using both first-principles calculations and the cluster expansion method. The cluster expansion method is used to describe the disordered alloys. It is revealed that, for both wurtzite (WZ) and zincblende (ZB) phase BexZn1−xO alloys, the formation enthalpies of all structures are positive for the whole range of composition. This implies the occurrence of miscibility gap when BeO and ZnO form alloys. A good comparison between the density functional theory used and the effective cluster interaction fitted formation enthalpies validates the cluster expansion method in the calculation of the formation enthalpies. The phase diagram of BexZn1−xO has been derived based on the theoretical calculations. It turns out that the inclusion of phonon contributions into the cluster expansion Hamiltonian affects markedly the substituent solubility of Be- and Zn-rich alloys. When lattice vibrations are considered, the solubility limits of Be in WZ-ZnO and Zn in WZ-BeO at 2000 K increase from 5.9% to 12.7% and from 0.7% to 3.8%, respectively, while the solubility limit of Be in ZB-ZnO reduces from 5.7% to 0.4% and that of Zn in ZB-BeO increases from 1.3% to 32.4%. A phase transition of BexZn1−xO from wurtzite to zincblende is predicted to occur around 1000 K.

Comparing the catalytic activity of the water gas shift reaction on Cu(3 2 1) and Cu(1 1 1) surfaces: Step sites do not always enhance the overall reactivity

A Density Functional Theory based first-principles kinetic Monte Carlo (kMC) study of the water gas shift reaction on the stepped Cu(321) surface is presented. We use the recently developed graph-theoretical kMC approach coupled with cluster expansion Hamiltonians to model the coverage-dependent energy barriers for the different surface processes, including adsorption/desorption, diffusion and other elementary chemical reactions, totalling 36 elementary steps, which allow two possible competitive mechanisms: surface redox and associative COOH. All results are compared to a previous kMC study on Cu(111). Both mechanisms are observed for Cu(321) surface with different extensions, whereas the associative COOH one was the dominant for Cu(111). The present study shows that, in spite of encompassing lower activation energy barriers, stepped surfaces do not necessarily have an overall larger catalytic activity. Coverage effects and the significant contribution of some of the reverse processes are behind this behaviour.

First-principles calculation of sulfur–selenium segregation in ZnSe1−xSx: The role of lattice vibration

First-principles phase diagram calculations based on density functional theory within the generalized gradient approximation in combination with Monte Carlo (MC) techniques and cluster expansion were performed for the ZnSe1−xSx alloys. Formation energies were used as a basis for fitting cluster expansion Hamiltonians. All formation energies of ZnSe1−xSx alloys are positive, showing that ZnSe1−xSx alloys is a miscibility gap system and has a tendency to phase separation. For ZnSe1−xSx alloys the phase diagram computed with conventional cluster expansion shows a miscibility gap with consolute temperature Tc=327K. The contributions of lattice vibration reduce Tc to 281K (about 15%). We presented a MC study of the spatial distribution of S and Se in ZnSe0.5S0.5 alloys. We found that, the system becomes more homogeneous including lattice vibration at lower temperature. It is consistence with the calculation of phase diagram of ZnSe1−xSx alloys.

Non-empirical phase equilibria in the Cr–Mo system: A combination of first-principles calculations, cluster expansion and Monte Carlo simulations

In this work, the low temperature part of the Cr–Mo phase diagram is calculated by a combination of first-principles calculations, cluster expansion, and Monte Carlo simulations without the use of empirical parameters. In brief, the miscibility gap in the bcc phase is computed through lattice Monte Carlo simulations with the energetics of the Cr–Mo system represented by cluster expansion Hamiltonians, wherein the Hamiltonians are firstly constructed using first-principles calculations of formation energies for supercells of various bulk compositions. It is demonstrated that the onset (critical) temperature of the miscibility gap is preliminarily computed to be Tc = 1720 K, and is significantly decreased to 1150 K once excess vibrational contribution to the free energy is considered. The resulting solid-state miscibility gap of the Cr–Mo binary system, largely improved by the contributions of atomic relaxation and positive excess volume of mixing in the Cr–Mo bcc alloys, also favors the approaches that extrapolate the Gibbs energy of the bcc phase retrieved at high temperatures, as in the CALPHAD method.

Phase behavior and mechanical properties of Ni–W studied by first-principles calculations and ab initio based thermodynamics

The substitutional ordering phenomena of face-centered cubic Ni-rich Ni–W alloys are investigated by means of cluster-expansion Hamiltonians based on density functional theory calculations. We observe a strong tendency for ordering in Ni–W characterized by the formation of tungsten chains along the 〈100〉 direction. While previous studies only report a D1a-type compound at a stoichiometry of Ni4W with this 〈100〉-ordering tendency, we predict that within the solubility range of tungsten in nickel, structures with tungsten 〈100〉-chains are the predominant form of short-range order even at elevated temperatures. In particular, we find a Pt8Ti-structured compound (which also shows this 〈100〉-ordering) to be stable in Ni–W at low temperatures and low tungsten contents. Even at high temperatures, Ni–W solid solutions are found to exhibit strong remnant order. The implications of these findings are discussed from a metallurgical perspective with an emphasis on the stiffening effect of tungsten alloying on Ni-rich matrices. Using these insights on the short-range order, we then determine a parametrization of the composition-dependent stiffening cij‾(xW) and ∂cij‾/∂xW.

Phase stability analysis of the InAs/GaAs (001) wetting layer from first principles

The large atomic-size mismatch between In and Ga and the large lattice-mismatch strain between InAs and GaAs make the InAs/GaAs (001) growth interface a complex alloy system, the understanding of which can enhance control of device synthesis and nanostructure self-assembly. We present a detailed first-principles analysis of the full progression of surface reconstructions observed on the InAs/GaAs(001) wetting layer during early stages of In deposition. We use systematic techniques to identify the most likely surface reconstruction prototypes of the InAs wetting layer on GaAs(001) using density functional theory (DFT) and then develop several cluster expansion Hamiltonians in order to thoroughly explore surface alloy disorder due to species substitution of In, Ga, and As at the surface. We use these results to construct a first principles 0-K surface phase diagram of the InAs wetting layer on GaAs(001) and test the sensitivity of our predicted phase diagram to DFT approximations and convergence errors. We find two alloy configurations of the $(4\ifmmode\times\else\texttimes\fi{}3)$ structural prototype that are likely ground-state surface reconstructions, and simulated scanning tunneling micrographs (STM) of these reconstructions indicate that they can explain prominent features of experimentally obtained STM of the InAs/GaAs $(4\ifmmode\times\else\texttimes\fi{}3)$ surface.

Parallel kinetic Monte Carlo simulation framework incorporating accurate models of adsorbate lateral interactions

Ab initio kinetic Monte Carlo (KMC) simulations have been successfully applied for over two decades to elucidate the underlying physico-chemical phenomena on the surfaces of heterogeneous catalysts. These simulations necessitate detailed knowledge of the kinetics of elementary reactions constituting the reaction mechanism, and the energetics of the species participating in the chemistry. The information about the energetics is encoded in the formation energies of gas and surface-bound species, and the lateral interactions between adsorbates on the catalytic surface, which can be modeled at different levels of detail. The majority of previous works accounted for only pairwise-additive first nearest-neighbor interactions. More recently, cluster-expansion Hamiltonians incorporating long-range interactions and many-body terms have been used for detailed estimations of catalytic rate [C. Wu, D. J. Schmidt, C. Wolverton, and W. F. Schneider, J. Catal. 286, 88 (2012)]. In view of the increasing interest in accurate predictions of catalytic performance, there is a need for general-purpose KMC approaches incorporating detailed cluster expansion models for the adlayer energetics. We have addressed this need by building on the previously introduced graph-theoretical KMC framework, and we have developed Zacros, a FORTRAN2003 KMC package for simulating catalytic chemistries. To tackle the high computational cost in the presence of long-range interactions we introduce parallelization with OpenMP. We further benchmark our framework by simulating a KMC analogue of the NO oxidation system established by Schneider and co-workers [J. Catal. 286, 88 (2012)]. We show that taking into account only first nearest-neighbor interactions may lead to large errors in the prediction of the catalytic rate, whereas for accurate estimates thereof, one needs to include long-range terms in the cluster expansion.

Strain effects on ensemble populations in AuPd/Pd(100) surface alloys

The effects of applied strain on the arrangement of atoms in AuPd∕Pd(100) surface alloys are studied using Monte Carlo simulations and cluster expansion Hamiltonians. The strain effects are found to be significant, with heteronuclear (Au-Pd) interactions more strongly enhanced by biaxial compression than homonuclear (Pd-Pd) ones. In particular, compressive strain causes an increase in the population of Pd monomers and second nearest-neighbor pairs of Pd monomers, both of which have been identified previously as important ensembles for various catalytic reactions. We also discuss the origin of these effects using density functional theory calculations of the surface electronic structure of strained AuPd∕Pd(100). Our findings may suggest an additional means of employing strain to tune the catalytic properties of surface alloys.

Thermodynamics of the Zr-O system from first-principles calculations

We investigate the electronic and thermodynamic properties of Zr and its oxides from first principles to elucidate phase stability in the Zr-O system. Hexagonally close-packed Zr is unusual in its ability to dissolve very high concentrations of oxygen over its interstitial octahedral sites, forming a variety of ordered suboxides that undergo both first-order and second-order phase transitions upon heating. We perform a first-principles, statistical-mechanical analysis of finite temperature phase stability of ZrO${}_{x}$ using a cluster expansion Hamiltonian and Monte Carlo calculations. This analysis predicts the existence of 0-K ground-state oxygen orderings at composition ZrO${}_{1/6}$, ZrO${}_{2/9}$, ZrO${}_{1/3}$, ZrO${}_{4/9}$, and ZrO${}_{1/2}$ along with evidence of an infinite sequence of ground-state suboxide orderings at intermediate oxygen concentrations consisting of different stackings of empty, $\frac{1}{3}$-filled and $\frac{2}{3}$-filled two-dimensional oxygen layers. We also predict the stability of a previously uncharacterized Zr-monoxide phase, which we label ${\ensuremath{\delta}}^{\ensuremath{'}}$-ZrO due to its crystallographic relation to $\ensuremath{\delta}$-TiO. The ${\ensuremath{\delta}}^{\ensuremath{'}}$-ZrO structure is equivalent to the high-pressure $\ensuremath{\omega}$-Zr phase but has interstitial oxygen ordering. Finally, as part of the technical implementation of our statistical mechanical study, we introduce a new algorithm to parametrize the coefficients of a cluster expansion Hamiltonian and apply a $k$-space analysis to rigorously track order-disorder phenomena at finite temperature.

Adsorbate-induced segregation: First-principles study for C/Pt25Rh75(100)

We develop a cluster-expansion (CE) formalism for the first-principles study of substrate-adsorbate interaction in order to study adsorbate-induced segregation. The simulation lattice, the configuration space, and the basis functions of the CE are bisected into separate but interacting subsystems. The method is applied to the C-contaminated Pt${}_{25}$Rh${}_{75}$$(100)$ surface. For such catalytic surfaces, it is essential to understand the change in the segregation profile of the substrate due to the presence of adsorbates. In accordance with experiments, we find that even a small amount of C impurities leads to a considerable decrease in the Pt segregation to the topmost substrate layer. We discuss the coupling of the substrate with the adsorbate layers in the CE Hamiltonian, the $T=0\phantom{\rule{4pt}{0ex}}\mathrm{K}$ stability diagram of C/Pt${}_{25}$Rh${}_{75}$$(100)$, and our prediction of the temperature-dependent segregation profile for different C contaminations.

Ab initiothermodynamics of intrinsic oxygen vacancies in ceria

Nonstoichiometric ceria(CeO${}_{2\ensuremath{-}\ensuremath{\delta}}$) is a candidate reaction medium to facilitate two-step water splitting cycles and generate hydrogen. Improving upon its thermodynamic suitability through doping requires an understanding of its vacancy thermodynamics. Using density functional theory (DFT) calculations and cluster expansion-based Monte Carlo simulations, we have studied the high-temperature thermodynamics of intrinsic oxygen vacancies in ceria. The DFT$+U$ approach was used to get the ground state energies of various vacancy configurations in ceria, which were subsequently fit to a cluster expansion Hamiltonian to efficiently model the configurational dependence of energy. The effect of lattice vibrations was incorporated through a temperature-dependent cluster expansion. Lattice Monte Carlo simulations using the cluster expansion Hamiltonian were able to detect the miscibility gap in the phase diagram of ceria. The inclusion of vibrational and electronic entropy effects made the agreement with experiments quantitative. The deviation from an ideal solution model was quantified by calculating as a function of nonstoichiometry, (a) the solid state entropy from Monte Carlo simulations, and (b) Warren-Cowley short range order parameters of various pair clusters.