Long-timescale simulations of diffusion in molecular solids

Kinetic processes play a crucial role in the formation and evolution of molecular layers. In this perspective we argue that adaptive kinetic Monte Carlo is a powerful simulation technique for determining key kinetic processes in molecular solids. The applicability of the method is demonstrated by simulating the diffusion of a CO admolecule on a water ice surface, which is an important process for the formation of organic compounds on interstellar dust grains. CO diffusion is found to follow Arrhenius behavior and the corresponding effective activation energy for diffusion is determined to be 50 ± 1 meV. A coarse graining algorithm is applied which greatly enhances the efficiency of the simulations at low temperatures, down to 10 K, without altering the underlying physical processes. Eventually, we argue that a combination of both on- and off-lattice kinetic Monte Carlo techniques is a good way for simulating large-scale processes in molecular solids over long time spans.     

Monte Carlo simulation of electrodeposition of copper: a multistep free energy calculation

Electrodeposition of copper (Cu) involves length scales of a micrometer or even less. Several theoretical techniques such as continuum Monte Carlo, kinetic Monte Carlo (KMC), and molecular dynamics have been used for simulating this problem. However the multiphenomena characteristics of the problem pose a challenge for an efficient simulation algorithm. Traditional KMC methods are slow, especially when modeling surface diffusion with large number of particles and frequent particle jumps. Parameter estimation involving thousands of KMC runs is very time-consuming. Thus a less time-consuming and novel multistep continuum Monte Carlo simulation is carried out to evaluate the step wise free energy change in the process of electrochemical copper deposition. The procedure involves separate Monte Carlo codes employing different random number criterion (using hydrated radii, bare radii, hydration number of the species, redox potentials, etc.) to obtain the number of species (CuCl(2) or CuSO(4) or Cu as the case may be) and in turn the free energy. The effect of concentration of electrolyte, influence of electric field and presence of chloride ions on the free energy change for the processes is studied. The rate determining step for the process of electrodeposition of copper from CuCl(2) and CuSO(4) is also determined.     

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.     

Spatially adaptive lattice coarse-grained Monte Carlo simulations for diffusion of interacting molecules

While lattice kinetic Monte Carlo (KMC) methods provide insight into numerous complex physical systems governed by interatomic interactions, they are limited to relatively short length and time scales. Recently introduced coarse-grained Monte Carlo (CGMC) simulations can reach much larger length and time scales at considerably lower computational cost. In this paper we extend the CGMC methods to spatially adaptive meshes for the case of surface diffusion (canonical ensemble). We introduce a systematic methodology to derive the transition probabilities for the coarse-grained diffusion process that ensure the correct dynamics and noise, give the correct continuum mesoscopic equations, and satisfy detailed balance. Substantial savings in CPU time are demonstrated compared to microscopic KMC while retaining high accuracy.     

Comparison of kinetic Monte Carlo and molecular dynamics simulations of diffusion in a model glass former

We present results from kinetic Monte Carlo (KMC) simulations of diffusion in a model glass former. We find that the diffusion constants obtained from KMC simulations have Arrhenius temperature dependence, while the correct behavior, obtained from molecular dynamics simulations, can be super-Arrhenius. We conclude that the discrepancy is due to undersampling of higher-lying local minima in the KMC runs. We suggest that the relevant connectivity of minima on the potential energy surface is proportional to the energy density of the local minima, which determines the "inherent structure entropy." The changing connectivity with potential energy may produce a correlation between dynamics and thermodynamics.     

Database of atomistic reaction mechanisms with application to kinetic Monte Carlo

Kinetic Monte Carlo is a method used to model the state-to-state kinetics of atomic systems when all reaction mechanisms and rates are known a priori. Adaptive versions of this algorithm use saddle searches from each visited state so that unexpected and complex reaction mechanisms can also be included. Here, we describe how calculated reaction mechanisms can be stored concisely in a kinetic database and subsequently reused to reduce the computational cost of such simulations. As all accessible reaction mechanisms available in a system are contained in the database, the cost of the adaptive algorithm is reduced towards that of standard kinetic Monte Carlo.     

Annealing contour Monte Carlo algorithm for structure optimization in an off-lattice protein model

We present a space annealing version for a contour Monte Carlo algorithm and show that it can be applied successfully to finding the ground states for an off-lattice protein model. The comparison shows that the algorithm has made a significant improvement over the pruned-enriched-Rosenbluth method and the Metropolis Monte Carlo method in finding the ground states for AB models. For all sequences, the algorithm has renewed the putative ground energy values in the two-dimensional AB model and set the putative ground energy values in the three-dimensional AB model.     

Halide adsorption on single-crystal silver substrates: dynamic simulations and ab initio density functional theory

We investigate the static and dynamic behaviors of a Br adlayer electrochemically deposited onto single-crystal Ag(100) using an off-lattice model of the adlayer. Unlike previous studies using a lattice-gas model, the off-lattice model allows adparticles to be located at any position within a two-dimensional approximation to the substrate. Interactions with the substrate are approximated by a corrugation potential. Using density functional theory (DFT) to calculate surface binding energies, a sinusoidal approximation to the corrugation potential is constructed. A variety of techniques, including Monte Carlo and Langevin simulations, are used to study the behavior of the adlayer. The lateral root-mean-square (rms) deviation of the adparticles from the binding sites is presented along with equilibrium coverage isotherms, and the thermally activated Arrhenius barrier-hopping model used in previous dynamic Monte Carlo simulations is tested.     

The XPK Package: A Comparison between the Extended Phenomenological Kinetic (XPK) Method and the Conventional Kinetic Monte Carlo (KMC) Method

Recently, we proposed the extended phenomenological kinetics (XPK) method, which overcomes the notorious timescale separation difficulty between fast diffusion and slow chemical reactions in conventional kinetic Monte Carlo (KMC) simulations. In the present work, we make a comprehensive comparison, based on the newly developed XPK package, between the XPK method and the conventional KMC method using a model hydrogenation reaction system. Two potential energy surfaces with different lateral interactions have been designed to illustrate the advantages of the XPK method in computational costs, parallel efficiency and the convergence behaviors to steady states. The XPK method is shown to be efficient and accurate, holding the great promise for theoretical modelling in heterogeneous catalysis, in particular, when the role of the lateral interactions among adsorbates is crucial.     

Achieving optimum hydrogen permeability in PdAg and PdAu alloys

The present work investigates both the diffusivity and permeability of hydrogen (H) in palladium-silver (PdAg) and palladium-gold (PdAu) alloys over a 400-1200 K temperature range for Pd(100-X)M(X), M=Ag or Au and X=0%-48% using density functional theory (DFT) and kinetic Monte Carlo simulations (KMC). DFT has been employed to obtain octahedral (O)-, tetrahedral (T)-, and transition state (TS)- site energetics as a function of local alloy composition for several PdAg and PdAu alloys with compositions in supercells of X=14.18%, 25.93%, 37.07%, and 48.15% with the nearest (NNs) and next nearest neighbors (NNNs) varied over the entire range of compositions. The estimates were then used to obtain a model relating the O, T, and TS energies of a given site with NN(X), NNN(X), and the lattice constant. The first passage approach combined with KMC simulations was used for the H diffusion coefficient predictions. It was found that the diffusion coefficient of H in PdAg alloy decreases with increasing Ag and increases with increasing temperature, matching closely with the experimental results reported in the literature. The calculated permeabilities of H in these novel binary alloys obtained from both diffusivity and solubility predictions were found to have a maximum at approximately 20% Ag and approximately 12% Au, which agree well with experimental predictions. Specifically, the permeability of H in PdAg alloy with approximately 20% Ag at 456 K is three to four times that of pure Pd, while the PdAu alloy at 12% Au is four to five times that of pure Pd at 456 K.     

Hybrid Monte Carlo implementation of the Fourier path integral algorithm

This paper formulates a hybrid Monte Carlo implementation of the Fourier path integral (FPI-HMC) approach with partial averaging. Such a hybrid Monte Carlo approach allows one to generate collective moves through configuration space using molecular dynamics while retaining the computational advantages associated with the Fourier path integral Monte Carlo method. In comparison with the earlier Metropolis Monte Carlo implementations of the FPI algorithm, the present HMC method is shown to be significantly more efficient for quantum Lennard-Jones solids and suggests that such algorithms may prove useful for efficient simulations of a range of atomic and molecular systems.     

Coarse-grained Simulations of Chemical Oscillation in Lattice Brusselator System

The oscillation behavior of a two-dimension lattice-gas Brusselator model was investigated. We have adopted a coarse-grained kinetic Monte Carlo (CG-KMC) procedure, where m￡m microscopic lattice sites are grouped together to form a CG cell, upon which CG processes take place with well-defined CG rates. Such a CG approach almost fails if the CG rates are obtained by a simple local mean field (s-LMF) approximation, due to the ignorance of correlation among adjcent cells resulting from the trimolecular reaction in this nonlinear system. By proper incorporating such boundary effects, thus introduce the so-called b-LMF CG approach. Extensive numerical simulations demonstrate that the b-LMF method can reproduce the oscillation behavior of the system quite well, given that the diffusion constant is not too small. In addition, the deviation from the KMC results reaches a nearly zero minimum level at an intermediate cell size, which lies in between the effective diffusion length and the minimal size required to sustain a well-defined temporal oscillation.     

A lattice Monte Carlo simulation of the FePt alloy using a first-principles renormalized four-body interaction

Although a lattice Monte Carlo method provides an effective, simple, and fast way to study thermodynamic properties of substitutional alloys, it cannot treat by itself the off-lattice effects, such as thermal vibrations and local distortions. Therefore, even if the interaction among atoms at lattice points is calculated accurately by means of first-principles calculations, the lattice Monte Carlo simulation overestimates the order-disorder phase transition temperature. In this paper, we treat this problem in the investigation of the FePt alloy, which has recently attracted considerable interest in its magnetic properties. We apply a simple version of the potential renormalization theory to determine the interaction among atoms, including partly the off-lattice effects by means of first-principles calculations. Then, we use the interaction to perform a lattice Monte Carlo simulation of the FePt alloy on a fcc lattice. From the results, we find that the transition temperature obtained after the present renormalization procedure becomes closer to the experimental value.     

Monte Carlo simulation of freely jointed chains with variable bond lengths

Monte Carlo simulation of freely jointed off-lattice chains with variable bond length is usually done with local random displacements of beads and with reptation moves (displacements of a bead along a chain). In dense systems, the acceptance ratio of reptations decreases strongly with density. We discuss versions of reptation moves, which are effective in dense systems. The idea, which comes from lattice systems, is to use a pseudovacancy (walker), which has the same size as a bead of a chain. The walker is attached to a neighbor chain and then another bead of that chain is cleaved. This is equivalent to a reptation move and a nonlocal displacement of the walker and since no free volume is needed, the move can be used with advantage in dense systems. A related technique are cooperative motions, which were introduced by T. Pakula for lattice models, where several chains change their conformation concomitantly. Such cooperative loops are implemented in the Monte Carlo algorithm by creating a temporary walker by cleaving a bead from a chain, moving it with reptations and finally annihilating the walker by attaching it to the same chain it was cleaved from. These moves and the condition of detailed balance are discussed in detail. As an example, we study the integrated autocorrelation time τ_int for the radius of gyration for a two-dimensional system. For reduced densities larger than 0,4, we find that with standard reptations and local bead displacements τ_int increases strongly with density. If reptations with either a permanent or a temporary walker are used in addition to local moves, the integrated autocorrelation time changes only very little with density and very dense systems can still be simulated efficiently.     

The influence of defects on the morphology of Si (111) etched in NH4F

We have implemented a kinetic Monte Carlo (KMC) simulation to study the effects of wafer miscut and wafer defects on the morphologies of Si (111) surfaces etched in NH4F. Although a conventional KMC simulation reproduced previously published results, it failed to produce the morphologies observed in our experiments. By introducing both dopant sites and lattice defect sites into the model, we are able to simulate samples having different dopant elements and densities as well as different defect concentrations. Using the modified KMC simulation, the simulated surface morphologies agree well with the morphologies observed in our experiments. The enhanced model also gives insights to the formation mechanism for multiple level stacking pits, a notable morphology on the etched surfaces of samples with very small miscut angles.     

Kinetic Monte Carlo study of submonolayer heteroepitaxial growth comparing Cu/Ni and Pt/Ni on Ni(100)

The surface patterns formed during submonolayer Cu/Ni and Pt/Ni heteroepitaxy upon a Ni(100) substrate have been investigated by kinetic Monte Carlo (KMC) simulations. The two-dimensional (2D) KMC simulations are based upon rate constants for a complete nearest-neighbor set of 729 uncorrelated Cu or Pt atoms and/or Ni site-to-site hopping mobilities. The rate constant activation energies are determined by classical-potential total-energy calculations using an embedded-atom method potential function from the literature. We find that diffusion of Cu atoms occurs at a faster rate and Pt atoms at a slower rate than that of Ni atoms on the flat Ni(100) surface, and the initial nucleation and growth patterns of 2D islands and the kinetic versus thermodynamic control of the growth vary as a consequence. In the temperature and deposition time regime in which we work, the Cu/Ni systems show less than random mixing, while the Pt/Ni systems show more than random mixing. The Cu/Ni system has bonding energies that result in a tendency to segregate toward subdomains of pure Ni and Cu, though kinetic effects in the epitaxy trap the development of the system at small subdomain sizes. The Pt/Ni system has bonding energies giving a tendency to intermix completely, while epitaxial kinetic effects modestly interfere with the complete mixing. The kinetically determined island morphologies under various Cu/Ni and Pt/Ni compositions and deposition rates differ substantially over time periods that are long on the deposition time scale, and therefore the island patterns can become frozen in place.     

Order parameters and molecular shapes of phosphonic acid dialkyl esters in a uniaxial hydrophobic environment. A theoretical study

The conformational behaviour of phosphonic acid dialkyl ester (PAE) molecules in a uniaxial hydrophobic environment was studied using the Monte Carlo method with the well-known Metropolis algorithm. The influence of the surrounding molecules on ordering processes of the intrinsic PAE molecules was taken into account by a molecular field. The intramolecular energy was calculated using 6–12 atom-atom and torsion potential functions. Conformations were analysed using torsion angle distributions, segment and bond order parameters. The order parameters of the C-H bonds are compared with experimental results of²H NMR. Further a general method of determining the effective shape of a molecule is presented. The shape is defined by probability distribution for finding atomic coordinates within volume elements during a Monte Carlo run. Thus the asymmetry of a molecule can be visualized. PAE molecules studied show a positive asymmetry in all cases.     

A graph-theoretical kinetic Monte Carlo framework for on-lattice chemical kinetics

Existing kinetic Monte Carlo (KMC) frameworks for the simulation of adsorption, desorption, diffusion, and reaction on a lattice often assume that each participating species occupies a single site and represent elementary events involving a maximum of two sites. However, these assumptions may be inadequate, especially in the case of complex chemistries, involving multidentate species or complex coverage and neighboring patterns between several lattice sites. We have developed a novel approach that employs graph-theoretical ideas to overcome these challenges and treat easily complex chemistries. As a benchmark, the Ziff-Gulari-Barshad system is simulated and comparisons of the computational times of the graph-theoretical KMC and a simpler KMC approach are made. Further, to demonstrate the capabilities of our framework, the water-gas shift chemistry on Pt(111) is simulated.     

Long time dynamics of a polymer with rigid body monomer units relating to a protein model: Comparison with the rouse model

Simulation data from an off-lattice polymer model are compared with data from the Rouse model. The polymer model is built of sequentially connected rigid monomer units that represent the amide planes of a protein backbone. The time propagation of the dynamics of the polymer model is performed by a Monte Carlo method. The elementary Monte Carlo steps correspond to local confomational changes in a window of three consecutive monomer units. The time autocorrelation functions of end-to-end distances from segments within the linear chain molecule are considered in detail. Both models exhibit a stretched exponential decay pattern. A comparison of the data from the Rouse model and the computer simulation provide an estimate of the time unit of 15 ps for a full scan of the algorithm for local conformational changes along the chain. With a conservative estimate of the parameters governing the Rouse model this time unit is four orders of magnitude longer than the elementary time step of a conventional computer simulation of polymer dynamics based upon the classical equations of motion for all atoms. © 1992 by John Wiley & Sons, Inc.     

Acceleration of Monte Carlo simulations through spatial updating in the grand canonical ensemble

A new grand canonical Monte Carlo algorithm for continuum fluid models is proposed. The method is based on a generalization of sequential Monte Carlo algorithms for lattice gas systems. The elementary moves, particle insertions and removals, are constructed by analogy with those of a lattice gas. The updating is implemented by selecting points in space (spatial updating) either at random or in a definitive order (sequential). The type of move, insertion or removal, is deduced based on the local environment of the selected points. Results on two-dimensional square-well fluids indicate that the sequential version of the proposed algorithm converges faster than standard grand canonical algorithms for continuum fluids. Due to the nature of the updating, additional reduction of simulation time may be achieved by parallel implementation through domain decomposition.