Kinetic Monte Carlo Scheme Research Articles

A combined kinetic Monte Carlo and phase field approach to model thermally activated dislocation motion

In this work, a kinetic Monte Carlo (KMC) scheme is integrated into a phase field dislocation dynamics (PFDD) model to account for thermal activation. While the implementation is general and could be applied to many different thermally activated processes, we apply the model to simulate a dislocation transmitting through an interface. Two cases are analyzed: a dislocation transmitting through a strong interface, and a dislocation transmitting through a more realistic grain boundary in tungsten. Results and trends for both cases are discussed in detail, and ultimately show that the combined KMC-PFDD algorithm is effective in capturing thermally activated dislocation processes.

Synthesis of gold nanoparticles in emulsion–emulsion precipitation route: Experiment, mechanism and simulation

Soft templates like water-in-oil emulsion drops (of 1-100μm diameter), have been used for synthesis of nanoparticles; with specific mean particle diameter and narrow particle size distribution. In this context, two separate emulsions have been reacted by an emulsion–emulsion precipitation route. Consequently, two different kinds of interacting populations coexist — emulsion drops evolving through drop coalescence and breakage, with multiple nanoparticle populations growing inside each emulsion drop. However, no model framework currently exists to simulate these two highly interacting populations, for predicting evolution of nanoparticle size. To this end, presently gold nanoparticles (GNPs) have been used as a model system, to explore the nanoparticle formation mechanism. The metal salt precursor and reducing agent (NaBH4 or N2H4) are prepared in two separate water-in-oil emulsions and are reacted upon mixing to form GNPs. A kinetic Monte Carlo (kMC) scheme accounting for drop coalescence, drop breakage, nucleation, particle growth and particle–particle coagulation was developed to handle simultaneous size and number evolution of aqueous emulsion drops and GNP populations. We find that, the mean GNP diameter decreases (from 6.9 to 3.5 nm) with increasing molar ratio (M) of reducing agent to precursor concentration. Further, increasing the precursor concentration at a fixed M, does not affect the final particle size. Our kMC simulation results not only agree well with experimental mean particle diameter and coefficient of variation (COV) of particle diameter, but was also validated with full particle size distribution of GNPs. So, M can be used as a tunable parameter for controlling GNP diameter experimentally.

Increase in rod diffusivity emerges even in Markovian nature.

Rod-shaped particles embedded in certain matrices have been reported to exhibit an increase in their center of mass diffusivity upon increasing the matrix density. This increase has been considered to be caused by a kinetic constraint in analogy with tube models. We investigate a mobile rodlike particle in a sea of immobile point obstacles using a kinetic Monte Carlo scheme equipped with a Markovian process, that generates gaslike collision statistics, so that such kinetic constraints do essentially not exist. Even in such a system, provided the particle's aspect ratio exceeds a threshold value of about 24, the unusual increase in the rod diffusivity emerges. This result implies that the kinetic constraint is not a necessary condition for the increase in the diffusivity.

On the microscopic behaviour of the vapour-liquid interface of methane-xenon mixture

Although there are many simulation studies of mixtures, there remains a gap in the understanding of the phase transitions occurring when molecules are progressively added to a closed system. When molecules are added, the fluid passes from a rarefied state through metastable states where liquid droplets are dispersed in a gas phase. The droplets increase in size and when a sufficient number of molecules is present in the system, the droplets coalesce and the system splits into vapour and liquid phases separated by a planar interface. Here we have explored the microscopic details of these phase transitions for a methane/xenon mixture at 189 K by performing canonical ensemble simulations using the kinetic Monte Carlo scheme for uniform and non-uniform systems. When vapour and liquid are in equilibrium, the interface separating the vapour and liquid phases is planar. We have investigated variations in pressure and composition, in the interfacial region, when mixtures of molecules of constant compositions are added into the system. These changes are different from those in single-component systems because the increase in the number of particles affects the pressure of the system and the compositions of the liquid and gas phases.

Nucleation of water clusters on functionalised graphite with kinetic Monte Carlo scheme

ABSTRACT Although activated carbon is generally hydrophobic, water can adsorb on its surface due to the presence of functional groups. During the initial interactions between functional groups and water molecules, atoms of the functional groups rotate around the sigma bonds to maximise the electrostatic interactions. This physical process was simulated within the framework of kinetic Monte Carlo (kMC). We extended kMC to account for the rotation of the atoms of functional groups around the sigma bonds by entropic sampling of the orientation space, and tested the modified kMC scheme for adsorption of water on graphite functionalised with either hydroxyl groups or carboxyl groups. Results show that the amount adsorbed is greater when the rotation of atoms around the sigma bonds is allowed. The interspacing of functional groups also affects the adsorbed density. It was found that the optimum distances between hydroxyl groups and carboxyl groups to initiate water cluster are 0.4 and 0.7 nm, respectively, because of the synergistic effect of electrostatic interactions between water molecule and adjacent functional groups. This effect is reduced when the spatial distance is either decreased (because of the increasing importance of the repulsion) or increased (because of the loss of the synergetic effect).

MOCOKI: A Monte Carlo approach for optimal control in the force of a linear kinetic model

A Monte Carlo framework for solving optimal control problems governed by kinetic models is presented. The focus is on a kinetic model with Keilson-Storer linear collision term and the control mechanism is an external space-dependent force. The purpose of this control is to drive an ensemble of particles to acquire a desired mean velocity and position and to reach a desired final configuration in phase space. For this purpose, a gradient-based computational strategy in the framework of Monte Carlo methods is developed. Results of numerical experiments successfully validate the proposed control framework. Program summaryProgram title: MOCOKICPC Library link to program files:https://doi.org/10.17632/6cb8fggwm2.1Licensing provisions: GNU General Public License 3Programming language: C/C++/PythonNature of problem: In many applications involving gases or plasma, it is not possible to assume a continuum and therefore models and methods that work at the mesoscopic level are needed. In this framework, kinetic models and Monte Carlo techniques play a central role. On the other hand, many present and envisioned application require to design control strategies by external forces that drive the evolution of these systems. For this purpose, it is necessary to implement optimal control techniques consistent with kinetic structures and Monte Carlo schemes. The novel methodology presented in this paper, although focusing on a linear kinetic model and a low-dimensional setting, can be extended in principle to nonlinear models and high-dimensional problems.Solution method: The proposed computational framework solves ensemble optimal control problems governed by linear kinetic models with a control-in-the-force mechanism. This framework combines transport techniques and Monte Carlo schemes for collision with numerical optimization methods.

Consistency of NVT, NPT, µVT and Gibbs (NV2T and NPT) with kinetic Monte Carlo schemes

We investigate the consistency of the simulation methodology for fluid phase equilibria and mixture adsorption from various ensembles [NVT (canonical), NPT (isothermal-isobaric systems), µVT (grand canonical) and Gibbs NV2T and Gibbs NPT] with the kinetic Monte Carlo (kMC) scheme. kMC is an alternative to Metropolis Monte Carlo (MC), in which chemical potential can be determined accurately. This is a major advantage in dense phases where the Widom method, usually employed in Metropolis MC, fails to provide sufficient accuracy. We have studied binary mixtures of methane with ethane to illustrate the consistency between the simulation methods for bulk fluids, and mixtures of ethanol–water and ammonia-water as examples of adsorption systems. The adsorption systems exhibit cooperative effects between water and ethanol and between water and ammonia: water by itself does not wet a carbon surface, but it adsorbs in large quantities when either ethanol or ammonia are present. Microscopic configurations taken from the simulations show a distinct difference between the adsorption mechanism in the water/ethanol mixtures and the water/ammonia mixtures.

A new interpretation of chemical potential in adsorption systems and the vapour–liquid interface

Chemical potential is a fundamental thermodynamic quantity that is constant everywhere in uniform or non-uniform systems at equilibrium. Because it is not a mechanical variable, its clear interpretation is elusive and its relationship to the energetics of the molecules that make up the system has not been established. In this work, we present a link between the chemical potential and molecular energetics, using a kinetic Monte Carlo scheme. We illustrate this new interpretation using argon as a model species giving examples for adsorption on a graphite surface and for a bulk vapour-liquid equilibrium (VLE). It was found that in either an adsorbed phase or a bulk liquid phase, the chemical potential is associated with repelling molecules, despite the number of these molecules being very small. In a rarefied phase it is associated with attracting molecules. In the interfacial regions in an adsorption system or in a VLE, the energetics of the repelling and attracting molecules contribute equally to the chemical potential.

A Multi-Scale Model for Insulin Self-Association Rates and Oligomerization Kinetics

Insulin, a therapeutic protein, can aggregate during manufacture and is considered detrimental to its quality. Early identification of aggregation inducing factors would be an essential aspect concerning the design of mitigation strategies during manufacturing processes. The time-scales and length-scales relevant for aggregation makes the use of all-atom (AA) explicit solvent simulations for prediction of aggregation thermodynamics and/or kinetics very difficult. Since aggregation is mainly driven by non-native interactions, we have used the physics based MARTINI coarse grained model to model protein-protein interactions. In an earlier work we had identified aggregation-prone, partially folded intermediates (PFIs) of insulin. Here we have employed a transient-complex model to determine self-association rate of native insulin (N) and the PFIs of insulin. To this end, a set of most probable conformations of the N-PFI complex were first obtained using docking and scoring on basis of MMPBSA energies. These docked complexes were then refined by separating the N and PFI monomers translationally followed by long molecular dynamics simulations. In two sets of test cases we show that the binding interface of the final N-PFI complex is similar in both the AA simulations and CG-MARTINI simulations with elastic network constraints. We hypothesize that this N-PFI complex can be taken as the diffusion limited transient state for the insulin homodimer on aggregation pathway. We use the TransComp server to calculate association rates for formation of this diffusion limited complex. The dimer association rates and monomer conformational transition rates are used to determine all the rates for formation of various small oligomeric species on insulin aggregation pathway. These rates are then input to a kinetic monte carlo (KMC) scheme to determine oligomerization kinetics of insulin.

Kinetic Monte Carlo model for 1-D migration in a field of strong traps: Application to self-interstitial clusters in W-Re alloys

In this work, we develop a kinetic Monte Carlo (KMC) scheme to describe 1-D migration in a field of traps. The model is parameterized solely on static calculations but is capable to predict the evolution of the system at finite temperature. The model is parameterized to describe the migration of self-interstitial clusters in the W-Re system that is of particular interest for fusion applications. The developed KMC model is benchmarked against molecular dynamics simulations and is found suitable to replace the latter at a fraction of the computational cost. We found that the presence of even traces of Re (starting from 0.1%) has a significant impact on the mobility of small 1-D migrating self-interstitial clusters. In addition to the computational scheme, the obtained results elucidate the impact of Re on void swelling and the suppression of void lattice formation under neutron irradiation of tungsten.

An efficient kinetic Monte Carlo scheme for computing Helmholtz free energy and entropy in bulk fluids and adsorption systems

We present an efficient kinetic Monte Carlo scheme to determine the Helmholtz free energy and entropy of bulk fluids and adsorption systems. The method is made possible because this technique enables the accurate determination of the chemical potential. The Helmholtz free energy, A, is obtained by integrating the chemical potential with respect to the number of molecules, at constant volume and temperature, and the entropy is then determined from the fundamental thermodynamic equation A = E − TS. The entropy of bulk argon is found to be in excellent agreement with values calculated from the established equation of state (EOS). In a system with two co-existing phases we show that our method can determine the surface tension at the Gibbs dividing surface, without recourse to the mechanical route of Irving, Kirkwood and Buff.We show that the intrinsic integral molecular Helmholtz free energy and entropy of the adsorbed phase, corrected for the surface excess, are independent of the size of the simulation box for all chemical potentials tested.The new procedure is illustrated for a range of gases commonly used in the characterization of porous solids (argon, nitrogen, carbon dioxide and ammonia) as adsorbates and a graphitic slit pore as the model adsorbent.

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.

Multiscale Computational Fluid Dynamics: Methodology and Application to PECVD of Thin Film Solar Cells

This work focuses on the development of a multiscale computational fluid dynamics (CFD) simulation framework with application to plasma-enhanced chemical vapor deposition of thin film solar cells. A macroscopic, CFD model is proposed which is capable of accurately reproducing plasma chemistry and transport phenomena within a 2D axisymmetric reactor geometry. Additionally, the complex interactions that take place on the surface of a-Si:H thin films are coupled with the CFD simulation using a novel kinetic Monte Carlo scheme which describes the thin film growth, leading to a multiscale CFD model. Due to the significant computational challenges imposed by this multiscale CFD model, a parallel computation strategy is presented which allows for reduced processing time via the discretization of both the gas-phase mesh and microscopic thin film growth processes. Finally, the multiscale CFD model has been applied to the PECVD process at industrially relevant operating conditions revealing non-uniformities greater than 20% in the growth rate of amorphous silicon films across the radius of the wafer.

A new kinetic Monte Carlo scheme with Gibbs ensemble to determine vapour–liquid equilibria

We present a new kinetic Monte Carlo scheme, as an alternative to the Gibbs ensemble Monte Carlo (GEMC) method, to determine vapour–liquid equilibria using a canonical ensemble in a system composed of two boxes. To illustrate the method, we have tested it with two systems: (1) argon over a range of temperatures from below the triple point to close to the critical point; (2) methane and ethane mixtures of various compositions at 180 K. The advantage of the new scheme is that chemical potentials of all components are accurately determined in both boxes. In particular, the chemical potential in the liquid box is determined much more accurately than with the Widom method employed in conventional GEMC simulations.

Monte-Carlo Simulation of the Ionic Transport of Electrolyte Solutions at High Concentrations Based on the Pseudo-Lattice Model

The theoretical understanding of the ionic transport in electrolyte solutions is not yet well-established at high concentrations, such as at C/M > 1. In our present study, two transport phenomena---self-diffusion and ionic conduction---of the electrolyte solution at high concentrations (roughly 0.05 ≤ C/M ≤ 3) are computationally simulated by a kinetic Monte-Carlo scheme. A swap mechanism in the three-dimensional pseudo-lattice is proposed to model the movement of ions and solvent, in which only the nearest-neighbor interaction is considered. The energy difference between the before- and after-swap states, based on which the stochastic Monte-Carlo process occurs, is found to be expressed in only two energetic terms; the coulombic repulsion between ions with the identical sign, +εii, and the heat of dissolution of the ionic crystal, Δdiss H. The self-diffusion coefficients of both ions and solvent decrease with the increase in the ionic concentration, which qualitatively agree with the experimental observation. The asymmetric bell-shape of the specific conductivity, σ, principally originates from the coulombic repulsion. The ionic concentration at which the maximum σ occurs well coincides with that of the experimentally observed results. The Δdiss H-dependency reveals that, ceteris paribus, σ marks the highest at around -2kT < Δdiss H< 0, out of the range of which σ is attenuated by either the ion-ion or ion-solvent interaction. Finally, the limitations in the model are addressed. The detail of the study is in press at the date of 15 April 2016 in J. Electrochem. Soc. Figure 1

Atomistic Simulations and Interfacial Morphology of Graphene Grown on SiC(0001) and SiC(000-1) Substrates

A Kinetic Monte Carlo scheme is applied to simulate with atomic resolution the synthesis of mono (few) layer(s) graphene (Gr) from a silicon carbide (SiC) substrate by selective evaporation of silicon (Si) atoms. The simulation computes the individual dynamics of the residual carbon (C) atoms which diffuse and reconfigure starting from the positions occupied in the SiC hexagonal lattice to the final Gr honeycomb structure. During the transition they gradually modify hybridization (from sp3 to sp2) and bond partners (from Si-C to C-C). We demonstrate that our method is able to recover the complex evolution steps of the epitaxial Gr on SiC in large systems for large time intervals. Moreover, the simulation results can be validated directly by means of comparison with experimental data when varying the material (e.g. initial surface configuration or polarity) or process (e.g. temperature and pressure) conditions.

Modeling p-type charge transport in thienoacene analogs of pentacene

The charge transport properties of two fused-ring thienoacenes, (a) the syn-isomer of dibenzo-thieno-dithiophene (DBTDT), packing in the solid state with a π–π stacking arrangement and also known as bis-benzo-thieno-thiophene (BBTT) and (b) C6-DBTDT, an alkylated derivative, packing in the more conventional herring-bone arrangement, are investigated computationally in the framework of the non-adiabatic hopping mechanism. Charge transfer rate constants are computed within the Marcus–Levich–Jortner formalism including a single effective mode treated quantum mechanically and are injected in a kinetic Monte Carlo scheme to propagate the charge carrier in the crystal. Charge mobilities are computed at room temperature with and without the influence of an electric field and are shown to compare very well with the measured mobilities in single-crystal devices. Both systems show an almost 1D charge transport with C6-DBTDT displaying about a ten times larger mobility value, in agreement with experiment. It is shown that the role of the HOMO-1 orbital is not relevant for BBTT, while it might contribute to a more marked 2D charge transport character for C6-DBTDT.

Monte-Carlo Simulation of the Ionic Transport of Electrolyte Solutions at High Concentrations Based on the Pseudo-Lattice Model

The theoretical understanding of the ionic transport in electrolyte solutions is not yet well-established at high concentrations, such as at C/mol · L− 1 > 1. In our present study, two transport phenomena—self-diffusion and ionic conduction—of the electrolyte solution at high concentrations (roughly 0.05 ⩽ C/mol · L− 1 ⩽ 3) are computationally simulated by a kinetic Monte-Carlo scheme. A “swap mechanism” in the three-dimensional pseudo-lattice is proposed to model the movement of ions and solvent, in which only the nearest-neighbor interaction is considered. The energy difference between the before- and after-swap states, based on which the stochastic Monte-Carlo process occurs, is found to be expressed in only two energetic terms; the coulombic repulsion between ions with the identical sign, , and the heat of dissolution of the ionic crystal, . The self-diffusion coefficients of both ions and solvent decrease with the increase in the ionic concentration, which qualitatively agree with the experimental observation. The asymmetric bell-shape of the specific conductivity, σ, principally originates from the coulombic repulsion. The ionic concentration at which the maximum σ occurs well coincides with that of the experimentally observed results. The -dependency reveals that, ceteris paribus, σ marks the highest at around , out of the range of which σ is attenuated by either the ion-ion or ion-solvent interaction. Finally, the limitations in the model are addressed.

An efficientmethod to determine chemical potential of mixtures in the isothermal and isobaric bulk phase with kineticMonte Carlo simulation

ABSTRACTA kinetic Monte Carlo (kMC) scheme in the NPT ensemble (constant number of molecules, pressure and temperature) has been developed to determine accurate chemical potentials for all components in a homogeneous mixture. The simulation requires two moves: (1) a displacement move and (2) a volume change move. In the former, the mobility rate of a selected molecule is determined by its interaction with all the other molecules in the system and is moved to a random position within the simulation box, according to the Rosenbluth algorithm, without any rejections (entropic sampling). The volume change move is decided by a comparison between either the instant pressure or the partial average pressure (with long-range correction) and the specified pressure and is carried out much less frequently than the displacement move. We applied this NPT scheme to a number of mixtures in both the gaseous and liquid phases, and show that the derived chemical potentials are accurate and reproducible. The method is recommended for obtaining chemical potentials in mixtures that are required as input in a grand canonical ensemble simulation.

Probing Potential Energy Surface Exploration Strategies for Complex Systems.

The efficiency of minimum-energy configuration searching algorithms is closely linked to the energy landscape structure of complex systems, yet these algorithms often include a number of steps of which the effect is not always clear. Decoupling these steps and their impacts can allow us to better understand both their role and the nature of complex energy landscape. Here, we consider a family of minimum-energy algorithms based, directly or indirectly, on the well-known Bell-Evans-Polanyi (BEP) principle. Comparing trajectories generated with BEP-based algorithms to kinetically correct off-lattice kinetic Monte Carlo schemes allow us to confirm that the BEP principle does not hold for complex systems since forward and reverse energy barriers are completely uncorrelated. As would be expected, following the lowest available energy barrier leads to rapid trapping. This is why BEP-based methods require also a direct handling of visited basins or barriers. Comparing the efficiency of these methods with a thermodynamical handling of low-energy barriers, we show that most of the efficiency of the BEP-like methods lie first and foremost in the basin management rather than in the BEP-like step.