When discussing algorithms used for modeling crystal growth, kinetic Monte Carlo emerges as a commonly employed method, such as its early application in atomic modeling of surface growth of GaAs using molecular beam epitaxy16. Additionally, Monte Carlo simulations of ‘units of growth’, which involve space-filling tiles or Voronoi polyhedra, present a coarse-grained alternative. This approach has proven applicable to a wide range of crystal types, encompassing porous crystalline materials, metal-organic frameworks, and ionic crystals17. On the other hand, in contexts such as high-pressure materials synthesis, where the focus lies solely on the final structure rather than the dynamical process of materials formation, commonly utilized approaches include genetic algorithms for molecular crystals18 and evolutionary algorithms for a broad range of materials systems19,20.In our study, we implement Monte Carlo simulations to examine the optimal structures formed from various diamondoid molecules. Constraints were imposed on the internal C-C bonds within each molecule, while those C-C bonds formed between the two molecules ‘were set to’ have the same bond angle and length as those within the bulk diamond lattice. Ab initio simulations predicted that the C-H bond dissociation energy was lower than that of the C-C bonds, so that dehydrogenation occurred prior to the C-C bond dissociation—an assumption shown reasonable for the formation of diamond using adamantane molecules2. Therefore, in our Monte Carlo simulations, only the carbon atoms and the C-C bond formation between diamondoid molecules are considered. Our algorithm can be regarded as a specific implementation of kinetic Monte Carlo, by taking advantage of the geometry of the diamondoids molecules and the isostructural carbon skeleton of diamondoids with nano-diamond. Here, the ‘unit of growth’ aligns with the carbon skeleton of the diamondoid molecules. Simultaneously, the problem is approached atomically, with each carbon occupying a site in a diamond lattice and the resulting C-C bonds being treated at the atomic level. Our algorithms leverage both the concept of the ‘unit of growth’ and atomic modeling methods, enabling efficient handling of large systems.The simulations were conducted within a closed box with dimensions of L × L × L (L = 40) in unit of diamond unit cell. To initiate the simulations, M molecules were randomly placed within the box. In order to expedite the initialization process, the i-th molecule was initialized to form at least one bond with the existing structures created by the preceding i-1 molecules. This initialization process ensured that the M molecules formed a connected structure.
Algorithm 1
Each Monte Carlo step in our molecular geometry specific simulation
V ← (xmax − xmin) × (ymax − ymin) × (zmax − zmin), where, for example, xmax is the maximum x coordinate of all the carbon atoms in the cluster.
\({{{{\mathcal{S}}}}}_{0}\leftarrow\) surface of the cluster. Surface consists of carbon atoms, which form <4 bonds with other carbon atoms.
Randomly pick a molecule M0 such that \({M}_{0}\cap {{{{\mathcal{S}}}}}_{0}\ne \varnothing\) and remove it.
n0 ← number of bonds broken when removing M0
\({{{{\mathcal{E}}}}}_{nn}\leftarrow\) empty nearest neighbor sites of the cluster after M0 is removed, i.e., for any atom that remains in the cluster, all its empty nearest neighbor sites belong to \({{{{\mathcal{E}}}}}_{nn}\)
Randomly pick an empty carbon atom site \({{{{\rm{C}}}}}_{1}\in {{{{\mathcal{E}}}}}_{nn}\). This guarantees the carbon atom with more surrounding empty sites has a higher probability to be selected to potentially form a C-C bond with atoms in M0.
counter ← 0
succeed ← False
(We will attempt “MAX_TRY = 10” times to put the molecule at the new position. “counter” records how many times we have tried. At each try, a random atom in M0 is chosen to be put at C1 and a random orientation of M0 around C1 is chosen. The possible orientations are restricted to those such that the C-C bonds of M0 align with those of the diamond lattice of the cluster).
repeat
Randomly pick an atom C0 ∈ M0.
Try to put C0 at the position of C1 and form bonds between C0 and the existing atoms in the cluster
if Bonds can form s.t. M0 does not overlap with existing atoms in the cluster then
succeed ← True
end if
counter ← counter + 1
until succeed ∣∣ counter ≥ MAX_TRY
if succeed then
n1 ← number of bonds formed after putting M0 to the new position
ΔNb ← n1 − n0
Calculate the new volume V and volume difference ΔV
p ← random number between 0 and 1
if \(p < \exp (-\beta (-{E}_{b}\Delta {N}_{b}+P\Delta V))\) then
M0 jump to the new position
end if
end if
if M0 does not jump to a new position then
put M0 back to the original position
end if
In each Monte Carlo step, we randomly select a diamondoid molecule to possibly update its position. The update of this diamondoid molecule position, or we call it a ‘movement’, breaks old C-C bonds and forms new C-C bonds with the remaining carbon atoms. (See Algorithm 1 for details) We calculated the enthalpy difference ΔH = -EbΔNb + PΔV before and after the movement of this diamondoid molecule, where ΔNb is the number of net carbon-carbon bonds formed between diamondoid molecules (if fewer C-C bonds formed than the C-C bonds breaking, ΔNb would be negative), Eb is the carbon-carbon bond formation energy, P represents the pressure, and V is the volume of the minimal cubic box that contains all the diamondoid molecules. To optimize the sampling efficiency of potential newly formed C-C bonds, we employ a weighted sampling approach for atoms in diamondoid structures as well as the subsequently formed connections. This molecular geometry specific approach ensures that C atoms surrounded by more empty sites are assigned higher weights during the sampling process. (See Algorithm 1 for details). The probability for which we accept the movement is \(\exp (-\beta \Delta H)\), where β = 1/kBT. When \(\exp (-\beta \Delta H)\ge 1\), we accept this movement; when \(0 < \exp (-\beta \Delta H) < 1\), we accept this movement with the probability \(\exp (-\beta \Delta H)\). The initial structure was further perturbed by setting a high initial temperature Ti, compared to the C-C bond dissociation energy. The temperature was then gradually reduced to a final low temperature Tf, at which we performed enough Monte Carlo simulations until the structure was converged. Specifically, 2 million Monte Carlo steps are performed at Ti, then temperature is linearly lowered to Tf in 2 million Monte Carlo steps. Finally, 10 million Monte Carlo steps are performed at Tf.
Molecular geometry specific Monte Carlo simulation of the efficacy of diamond crystal formation from diamondoids
