Crystal Growth Kinetic Monte Carlo (cgkmc)

• Crystal Growth Kinetic Monte Carlo (cgkmc) is a class of stochastic, event-driven simulations that model crystal growth by executing discrete elementary events to bridge microscopic kinetics and macroscopic morphology.
• It employs diverse surface representations—such as SOS height fields and graph-based models—and tailored event catalogs to simulate phenomena from protein crystallization to thin-film epitaxy.
• Recent extensions of cgkmc integrate off-lattice diffusion, coarse-grained formulations, and machine-learning augmentation, addressing challenges like rate uncertainties and transport coupling.

Searching arXiv for recent and foundational papers on crystal growth kinetic Monte Carlo and cgkmc. Searching for the 2010 roughening-transition paper and related crystal-growth KMC work. Crystal Growth Kinetic Monte Carlo (cgkmc) denotes a class of stochastic, event-driven simulation methods for crystal growth in which the evolving crystal surface is represented explicitly and advanced through discrete elementary events such as adsorption, desorption, diffusion, attachment, detachment, and nucleation. In the literature summarized here, the term spans both general surface-growth KMC formulations and a specific open-source Python package, cgkmc, for morphology prediction at solid–liquid interfaces (Lutsko et al., 2010, Jeffries et al., 29 Sep 2025). Across these uses, cgkmc connects microscopic energetics and local kinetics to macroscopic observables such as growth velocity, island density, roughness, facet selection, texture evolution, and growth-mode transitions. The method has been applied to protein-crystal surfaces, molecular crystals such as PETN, confined growth and dissolution, electrodeposition, thin-film epitaxy, two-dimensional materials, and coarse-grained PECVD film growth (Lutsko et al., 2010, Jeffries et al., 29 Sep 2025, Høgberget et al., 2020, Treeratanaphitak et al., 2013, Xu et al., 2016, Aldana et al., 2024, Wu et al., 2018, Budagosky et al., 2021).

1. Definition and conceptual scope

cgkmc models crystal growth as a continuous-time stochastic process on a discrete representation of the surface. The representation may be a solid-on-solid (SOS) height field, a lattice occupation model, a graph of atomistic lattice sites, or a coarse-grained three-dimensional lattice whose sites encode mesoscopic state variables such as orientation (Lutsko et al., 2010, Jeffries et al., 29 Sep 2025, Høgberget et al., 2020, Budagosky et al., 2021). The central abstraction is that the crystal morphology emerges from repeated execution of local events whose probabilities are proportional to kinetic rates.

The 2010 study of protein-crystal growth formulated a reduced SOS KMC model and showed that a transition between smooth step-flow and rough multi-nucleation growth can arise even when the two-dimensional nucleation barrier remains finite (Lutsko et al., 2010). The 2025 PETN work defined cgkmc more narrowly as a graph-based, atomistic framework for non-equilibrium morphology prediction at solid–liquid interfaces and implemented it in the open-source Python package cgkmc (Jeffries et al., 29 Sep 2025). Other papers treated related KMC formulations as cgkmc in practice even when they used different code bases or application-specific parameterizations, including off-lattice diffusion coupling in confinement, EAM-based electrodeposition, HCP-lattice GaN homoepitaxy, and substrate-sensitive growth of TMD monolayers (Høgberget et al., 2020, Treeratanaphitak et al., 2013, Xu et al., 2016, Aldana et al., 2024, Wu et al., 2018).

A common feature across these formulations is rejection-free residence-time dynamics. In the standard form, the total event rate is $R = \sum_i k_i$, one event is sampled with probability $k_i/R$, and the simulation clock advances by $\Delta t = -\ln(u)/R$, with $u \in (0,1)$ uniformly distributed (Lutsko et al., 2010, Jeffries et al., 29 Sep 2025, Høgberget et al., 2020, Aldana et al., 2024, Wu et al., 2018, Xu et al., 2016, 1904.02561). This direct mapping between event catalog and physical time is one reason cgkmc is used to bridge atomistic kinetics and experimentally accessible growth regimes.

2. Surface representations and elementary processes

The surface representation used in cgkmc depends on the physical system. The protein-crystal study used a one-dimensional SOS height field $h_i$ on a lattice with periodic boundary conditions, while noting that cgkmc is typically performed on 2D lattices representing terraces (Lutsko et al., 2010). The PETN morphology model represented the crystal supercell as a “colored,” undirected graph $G=(V,Q)$, with binary site occupations $x_i \in \{0,1\}$ and pairwise interactions stored in a symmetric matrix $Q$ (Jeffries et al., 29 Sep 2025). Confined growth and dissolution used a $(2+1)$-D SOS height field $h(x,y)$ on a cubic lattice coupled directly to off-lattice solute diffusion (Høgberget et al., 2020). GaN homoepitaxy employed an ideal HCP lattice of Ga sites, whereas the MoS$k_i/R$0 and WS$k_i/R$1 studies used on-lattice models on hexagonal or substrate-sensitive lattices to resolve orientation-dependent edge kinetics (Xu et al., 2016, Aldana et al., 2024, Wu et al., 2018).

The event catalog also varies with the growth problem. In general surface-growth cgkmc, the event set may include adsorption or attachment to terraces, desorption or detachment, surface diffusion, edge diffusion, step-edge attachment, step-edge detachment, two-dimensional island nucleation, and defect-mediated step generation (Lutsko et al., 2010). The PETN implementation restricted the catalog to adsorption from solution and evaporation at solid–liquid interfaces, omitting surface diffusion and solvation effects in the present work (Jeffries et al., 29 Sep 2025). The confined-growth model combined nearest-neighbor dissolution on the crystal surface with MFPT-based deposition events supplied by off-lattice free-solute diffusion (Høgberget et al., 2020). The Cu electrodeposition model included deposition plus three surface diffusion mechanisms—terrace hopping, atom exchange, and step-edge atom exchange—under galvanostatic conditions (Treeratanaphitak et al., 2013). In MoS$k_i/R$2, the event set explicitly resolved adsorption, desorption, terrace diffusion, edge attachment, and edge migration; nucleation and coalescence were not simulated because the calculation began from a single triangular seed (Aldana et al., 2024).

A compact comparison of representative cgkmc formulations is useful because the same term covers materially different implementations.

System	Surface representation	Principal events
Protein-crystal surfaces	1D SOS height field $k_i/R$3	Adsorption, evaporation, critical-cluster stabilization, defect-driven step sources (Lutsko et al., 2010)
PETN morphology prediction	Binary occupation vector on graph $k_i/R$4	Adsorption and evaporation at interfacial sites (Jeffries et al., 29 Sep 2025)
Confined growth/dissolution	$k_i/R$5-D SOS plus off-lattice solute diffusion	Dissolution and MFPT-based deposition (Høgberget et al., 2020)
Cu electrodeposition	On-lattice FCC $k_i/R$6 surface with EAM energetics	Deposition, hopping, atom exchange, step-edge exchange (Treeratanaphitak et al., 2013)
Monolayer MoS$k_i/R$7	Hexagonal lattice	Adsorption, desorption, hopping, edge attachment, edge migration (Aldana et al., 2024)

This diversity implies that cgkmc is best understood as a methodological family rather than a single canonical model. A plausible implication is that model transferability depends less on the label “cgkmc” than on whether the chosen state space and event set resolve the kinetic bottlenecks relevant to the target material system.

3. Rate laws, energetics, and time evolution

The most common rate law in cgkmc is Arrhenius. In the standard thermally activated form used in several of the cited studies,

$k_i/R$8

where $k_i/R$9 is an attempt frequency and $\Delta t = -\ln(u)/R$0 is the event-specific activation barrier (Lutsko et al., 2010, Aldana et al., 2024, Wu et al., 2018, 1904.02561). This form underlies terrace diffusion, edge diffusion, detachment, edge migration, and step crossing in many crystal-growth KMC models.

Different works parameterize $\Delta t = -\ln(u)/R$1 differently. In the protein-crystal formulation, activation barriers encode migration, bond breaking, attachment penalties, and nucleation energetics, with supersaturation $\Delta t = -\ln(u)/R$2 biasing attachment and detachment through detailed-balance-consistent barriers (Lutsko et al., 2010). In the PETN package, the total interaction energy is

$\Delta t = -\ln(u)/R$3

and flipping the occupation at site $\Delta t = -\ln(u)/R$4 gives local energy changes

$\Delta t = -\ln(u)/R$5

which are then inserted into the kinetic model (Jeffries et al., 29 Sep 2025). Evaporation is treated by harmonic transition-state theory as

$\Delta t = -\ln(u)/R$6

while adsorption from solution is modeled by diffusion to a spherical sink, yielding a site-independent interfacial deposition rate

$\Delta t = -\ln(u)/R$7

with $\Delta t = -\ln(u)/R$8 the instantaneous effective crystal radius (Jeffries et al., 29 Sep 2025).

In the confined-growth model, the deposition rate is not imposed phenomenologically but derived from a mean-first-passage-time estimate. For a solute particle $\Delta t = -\ln(u)/R$9 at distance $u \in (0,1)$0 from a target top site $u \in (0,1)$1, the rate is

$u \in (0,1)$2

with radial and pathfinding variants for the definition of $u \in (0,1)$3 (Høgberget et al., 2020). This paper showed that the $u \in (0,1)$4 scaling follows from a diffusive MFPT $u \in (0,1)$5 modulated by the small target solid-angle fraction $u \in (0,1)$6 (Høgberget et al., 2020).

Other studies use system-specific energetics. The Cu electrodeposition model computes the energy difference $u \in (0,1)$7 from an EAM Hamiltonian and sets the site propensity for hopping, exchange, or step-edge exchange to

$u \in (0,1)$8

and

$u \in (0,1)$9

with mechanism-specific barriers $h_i$0, $h_i$1, and $h_i$2 (Treeratanaphitak et al., 2013). The GaN HCP model uses a nearest-neighbor bond-counting site energy $h_i$3 and a jump rate $h_i$4, with the activation energy depending on whether the move is uphill or downhill in local energy (Xu et al., 2016).

The common rejection-free time integrator is the residence-time or BKL algorithm. The PETN code, the MoS$h_i$5 implementation, the WS$h_i$6 study, the confined-growth solver, and the GaN and general lattice-KMC treatments all use the same essential update: $h_i$7 after selecting an event with probability proportional to its rate (Jeffries et al., 29 Sep 2025, Aldana et al., 2024, Wu et al., 2018, Høgberget et al., 2020, Xu et al., 2016, 1904.02561). This formal commonality masks substantial physical differences in how rates are constructed.

4. Nucleation, step flow, and kinetic roughening

A central problem in cgkmc is the competition between step-flow growth and terrace nucleation. The 2010 protein-crystal study provides the clearest explicit formulation. For a circular terrace island of radius $h_i$8, the reversible work of formation is

$h_i$9

with critical radius and barrier

$G=(V,Q)$0

The paper emphasized that a roughening-like transition need not require $G=(V,Q)$1; a finite nucleation barrier can still coexist with a kinetic crossover to multi-nucleation growth when supersaturation increases the probability current over the barrier (Lutsko et al., 2010).

In that study, low supersaturation produced smooth surfaces whose growth proceeded by step flow driven by steps nucleated heterogeneously at defects such as dislocations or “wall-like” features. At high supersaturation, many islands nucleated homogeneously across terraces, the surface became rough, and the growth velocity increased more rapidly with supersaturation than in the smooth regime (Lutsko et al., 2010). The key control variable was the competition between the step-sweep time and the nucleation time. With terrace width $G=(V,Q)$2, step velocity $G=(V,Q)$3, nucleation rate $G=(V,Q)$4, and terrace area $G=(V,Q)$5,

$G=(V,Q)$6

and the dimensionless ratio

$G=(V,Q)$7

signals the crossover near $G=(V,Q)$8 (Lutsko et al., 2010).

The kinetic signature of this crossover is not unique to protein crystals. In submonolayer MBE, island density and morphology similarly depend on the competition between deposition and surface transport. For irreversible aggregation with critical nucleus size $G=(V,Q)$9, the simple-cubic MBE study used the scaling law

$x_i \in \{0,1\}$0

with $x_i \in \{0,1\}$1 for the isotropic case and $x_i \in \{0,1\}$2 for anisotropic substrate diffusion (Devkota et al., 2019). The GaN homoepitaxy study extracted island spacing laws of the form

$x_i \in \{0,1\}$3

with exponents $x_i \in \{0,1\}$4 for m-plane and $x_i \in \{0,1\}$5 for c-plane, leading to power-law island-spacing exponents of $x_i \in \{0,1\}$6 and $x_i \in \{0,1\}$7, respectively (Xu et al., 2016).

These results collectively show that cgkmc treats roughening not as a single universality class but as a kinetic outcome of competing timescales: arrival, diffusion, step propagation, edge relaxation, and nucleation. This suggests that “kinetic roughening” in cgkmc applications is often operationally identified through observables such as roughness, island density, diffuse scattering, or slope changes in growth velocity rather than through a single microscopic mechanism.

5. Morphology prediction, anisotropy, and materials-specific implementations

One major branch of cgkmc aims at direct morphology prediction. The 2025 PETN study is the clearest example. Its graph-based atomistic routine predicts crystal shape by evolving the discrete occupancy of a supercell using only solid-phase intermolecular interaction parameters and a diffusion-based adsorption model for the liquid phase (Jeffries et al., 29 Sep 2025). For PETN, the predicted dominant planes $x_i \in \{0,1\}$8 and $x_i \in \{0,1\}$9 and the evolving surface energy

$Q$0

agree closely with attachment-energy-based predictions, while a constant $Q$1 yields an unphysical rod-like crystal dominated by $Q$2 (Jeffries et al., 29 Sep 2025). In this formulation, morphology is not imposed by a Wulff construction; it emerges from the spatial competition between adsorption and evaporation rates.

A second major branch uses cgkmc to study anisotropy in 2D materials. For monolayer MoS$Q$3, a hexagonal-lattice model incorporated substrate effects through local changes in adsorption rate $Q$4, desorption barrier $Q$5, diffusion barrier $Q$6, and edge-migration barrier $Q$7 (Aldana et al., 2024). At $Q$8 with attempt frequency $Q$9, the study used orientation-dependent edge-attachment barriers

$(2+1)$0

so ZZ edges advance faster than AC edges and triangular shapes with ZZ facets emerge (Aldana et al., 2024). The work identified a practical edge-reshaping Damköhler-like number,

$(2+1)$1

with compact triangles favored for $(2+1)$2 and branched morphologies for $(2+1)$3 (Aldana et al., 2024).

The WS$(2+1)$4 study on ST-X quartz encoded substrate anisotropy through a binary adsorption landscape of strong “Si” and ordinary “O” domains, with DFT-derived adsorption-energy contrast $(2+1)$5 (Wu et al., 2018). It found that, besides variation of substrate adsorption, the chalcogen-to-metal ratio $(2+1)$6 is a major contributor to anisotropic growth and to the polarization of undergrowth and overgrowth. At $(2+1)$7, the simulations identified two linear $(2+1)$8-$(2+1)$9 relations,

$h(x,y)$0

for isotropy and

$h(x,y)$1

for high anisotropy (Wu et al., 2018). Shape anisotropy was quantified through the anisotropic growth ratio (AGR), with $h(x,y)$2 corresponding to isotropic growth and $h(x,y)$3 to extreme anisotropy (Wu et al., 2018).

For GaN homoepitaxy, morphology was linked to crystallographic orientation through effective step-edge energies and ES barriers on an HCP lattice (Xu et al., 2016). The simulations showed that elongated islands on m-plane surfaces arise primarily from anisotropy of step-edge energy rather than adatom diffusion anisotropy (Xu et al., 2016). For coarse-grained PECVD of TiO$h(x,y)$4 and ZnO, orientation and local structural descriptors were embedded directly into the coarse particles. In ZnO, each lattice particle carried a discretized orientation label, and the orientation-dependent local barrier term was built from an anisotropy ratio

$h(x,y)$5

together with an orientation-mismatch weakening factor

$h(x,y)$6

allowing simulation of texture competition between [002], [101], and [100] growth (Budagosky et al., 2021).

Across these systems, cgkmc serves not merely as a growth-rate calculator but as a morphology generator. The specific observables differ—dominant facets in PETN, compact versus branched flakes in MoS$h(x,y)$7, AGR in WS$h(x,y)$8, reciprocal-space diffuse satellites in GaN, texture fractions in ZnO—but each arises from local event asymmetries encoded into the rate model.

6. Transport coupling, multiscale extensions, and machine-learning augmentation

Several works extend cgkmc beyond conventional on-lattice surface kinetics. The confined-growth and dissolution model directly coupled an SOS surface KMC solver to a free-solute diffusion model in which liquid particles move off-lattice by Gaussian random walks with variance $h(x,y)$9 in each Cartesian direction (Høgberget et al., 2020). The coupling was achieved without a continuum diffusion solver: deposition rates were computed directly from MFPT estimates, and all nonparticipating solute particles were advanced by off-lattice diffusion during each KMC time step (Høgberget et al., 2020). The paper found that free-diffusion MFPT models reproduce the macroscopic limit of surface dynamics, perform well in tight confinement, and yield lower equilibrium roughness and shorter height autocorrelation lengths than uniform-concentration and lattice-diffusion reference models (Høgberget et al., 2020).

A different extension appears in coarse-grained PECVD growth, where the atomic-scale surface is replaced by coarse particles on a cubic lattice. There the activation energy of a surface hop is written as

$k_i/R$00

with the fluctuating term $k_i/R$01 computed on the fly from the local mesoscale morphological or structural landscape (Budagosky et al., 2021). For TiO$k_i/R$02, $k_i/R$03 was drawn from a discretized interval $k_i/R$04 and spatially correlated through neighbor averaging, while for ZnO it was orientation-dependent (Budagosky et al., 2021). This formulation allowed simulation of shadowing-driven columnar growth, roughening exponents, and texture evolution under PECVD conditions not accessible to atomistic KMC at comparable scales (Budagosky et al., 2021).

The most explicit atomistic augmentation of cgkmc with machine learning appears in the 2026 self-evolving ML-based method for Ag/Ag$k_i/R$05 thin-film growth (Kimari et al., 29 May 2026). That framework builds Gaussian-process regression models over SOAP descriptors of local atomic environments during runtime and triggers NEB calculations when the predictive uncertainty exceeds $k_i/R$06 (Kimari et al., 29 May 2026). Barriers are then predicted by ML once sufficient training data accumulate, so expensive NEB calculations are progressively replaced. The method uses a lattice-free KMC on a dynamically refined adsorption-site network, with deposition at $k_i/R$07 monolayers/s and a constant diffusion prefactor $k_i/R$08 (Kimari et al., 29 May 2026). It reproduces terrace, edge, corner, step-ascent, and step-descent barriers in close agreement with NEB and literature values, while predicting island densities and morphologies consistent with experimental and theoretical expectations (Kimari et al., 29 May 2026).

These extensions indicate that cgkmc now spans a broad spectrum of multiscale strategies: explicit diffusion coupling, coarse-grained mesoscale event models, graph-based atomistic morphology prediction, and uncertainty-aware ML acceleration. This suggests that the methodological frontier of cgkmc is not a single improved lattice Hamiltonian but the systematic coupling of stochastic surface kinetics to more faithful descriptions of transport, energetics, and local environment.

7. Assumptions, limitations, and methodological issues

Despite its breadth, cgkmc remains strongly model-dependent. A recurring assumption is the use of SOS or lattice constraints. The protein-crystal, confined-growth, MBE, and many lattice-KMC implementations disallow overhangs and often neglect bulk vacancies or long-range relaxation (Lutsko et al., 2010, Høgberget et al., 2020, Aryal, 2024, Devkota et al., 2019, 1904.02561). The Cu electrodeposition model is on-lattice and uses fixed activation barriers from the literature rather than locally computed saddle-point barriers, even though local energetic bias is evaluated with an EAM Hamiltonian (Treeratanaphitak et al., 2013). The PETN package neglects surface diffusion and solvent structure in the present implementation (Jeffries et al., 29 Sep 2025). The MoS$k_i/R$09 model neglects nucleation, coalescence, explicit Mo kinetics, and multilayer growth (Aldana et al., 2024). The coarse-grained PECVD model does not include detailed plasma-surface chemistry, ion-induced resputtering, or true hexagonal facets (Budagosky et al., 2021).

Rate uncertainty is another general issue. The practical review of lattice KMC emphasizes that missing elementary processes or errors in rate constants can alter kinetics by orders of magnitude because the rates depend exponentially on activation barriers (1904.02561). The electrodeposition study provides a concrete illustration of a compromise: realistic multi-body energetics from EAM are combined with fixed mechanism-specific barriers $k_i/R$10, $k_i/R$11, and $k_i/R$12 (Treeratanaphitak et al., 2013). The ML-augmented Ag/Ag$k_i/R$13 method addresses this issue by computing barriers on the fly with NEB when uncertainty is large, but it still assumes a constant attempt frequency across diffusion events and excludes concerted moves beyond the enumerated catalog (Kimari et al., 29 May 2026).

Transport and environmental dependence also remain difficult. The PETN study notes that under its target-size control scheme, $k_i/R$14 and $k_i/R$15 primarily affect timescale rather than final morphology, which conflicts with experimental solvent dependence and motivates more realistic prefactor parameterizations (Jeffries et al., 29 Sep 2025). The confined-growth paper shows that line-of-sight effects can often be ignored near equilibrium or during dissolution, but become increasingly important at higher growth speed (Høgberget et al., 2020). The MoS$k_i/R$16 study shows that substrate etching can be represented by local changes in $k_i/R$17, $k_i/R$18, and $k_i/R$19, yet this remains an effective description rather than an explicit treatment of surface chemistry (Aldana et al., 2024).

Finally, there is no single “cgkmc metric set.” Depending on the application, researchers analyze growth velocity $k_i/R$20, roughness $k_i/R$21, slope roughness $k_i/R$22, island density $k_i/R$23, island-size distributions, autocorrelation length $k_i/R$24, AGR, texture coefficients, reciprocal-space anti-Bragg intensities, or surface energies (Lutsko et al., 2010, Jeffries et al., 29 Sep 2025, Høgberget et al., 2020, Wu et al., 2018, Xu et al., 2016, Budagosky et al., 2021). This diversity is methodologically appropriate but can obscure cross-study comparison.

A common misconception is that cgkmc is synonymous with a single package or a single lattice-gas formalism. The literature does not support that view. The PETN package is one concrete software realization named cgkmc (Jeffries et al., 29 Sep 2025), but the broader concept encompasses SOS growth models, graph-based atomistic morphology prediction, EAM-driven electrodeposition, MFPT-coupled diffusion-growth solvers, and ML-augmented dynamic-site KMC (Lutsko et al., 2010, Treeratanaphitak et al., 2013, Høgberget et al., 2020, Kimari et al., 29 May 2026). Another misconception is that roughening transitions necessarily imply the vanishing of the two-dimensional nucleation barrier; the protein-crystal study explicitly demonstrates a roughening-like transition at finite nucleation barrier (Lutsko et al., 2010).

In that broader sense, cgkmc is best characterized as a rigorously time-resolved mesoscale simulation framework in which the fidelity of predictions depends on the physical adequacy of the event catalog, the barrier model, and the coupling to transport. The literature suggests that its principal strength lies in exposing how local stochastic kinetics produce emergent morphology, and its principal limitation lies in the difficulty of parameterizing those kinetics without either oversimplification or prohibitive computational cost.