Molecular geometry specific Monte Carlo simulation of the efficacy of diamond crystal formation from diamondoids (2024)

Diamondoids are a series of hydrogen-terminated nanometer-sized hydrocarbons where the carbon skeletons represent the nano-diamond structures¹. Various diamondoid molecules, including the lower diamondoids (adamantane C₁₀H₁₆, diamantane C₁₄H₂₀, and triamantane C₁₈H₂₄) and higher diamondoids (C₂₂H₂₈ and higher polymantanes), have been utilized to synthesize diamond crystals using high pressure high temperature (HPHT)^2,3,4,5,6 or chemical vapor deposition (CVD)^5,7,8,9 methods, in the quest of finding better materials to make higher quality diamond crystals at a reduced cost.

In the HPHT experiments, diamond from lower diamondoid molecules forms at much lower temperatures compared with conventional (hydro)carbon allotropes at a comparable pressure, suggesting a reduced transformation energy barrier, owing to the structural similarities between diamondoids and the bulk diamond². Formation of diamond crystals from diamondoids at HPHT usually involves hydrogen dissociation². The hydrogen-stripped carbon skeletons then reorient to form carbon-carbon bonding that builds into larger diamond crystals. The diamondoids-to-diamond formation can be seen as a three-dimensional origami process for the carbon skeletons. Since different diamondoid molecules have different geometries of the carbon skeletons, specific types of diamondoids can offer the advantage of constructing diamond crystals without supplemental breaking of internal carbon-carbon bonds; while other types of diamondoid molecules form lower quality diamond with increased carbon vacancies due to increased necessity to reorient the carbon skeletons.

In order to aidin the experimental search for the optimal candidates and conditions for diamond synthesis, we performed a Monte Carlo simulation to study the potentials of different diamondoids in constructing diamond crystals. Six diamondoid molecules are studied: adamantane C₁₀H₁₆, diamantane C₁₄H₂₀, triamantane C₁₈H₂₄, [121]tetramantane C₂₂H₂₈, [1(2,3)4]pentamantane C₂₆H₃₂, and [1212]pentamantane C₂₆H₃₂. Our Monte Carlo results show that higher diamondoids have a higher efficacy to form diamond crystals compared with lower diamondoids, given the carbon skeleton for each diamondoid molecule remains intact without breaking internal C-C bonds during materials synthesis.

Monte Carlo simulations were conducted on six diamondoid molecules, including adamantane, diamantane, triamantane, [121]tetramantane, [1(2,3)4]pentamantane, and [1212]pentamantane, with the final structure of diamond formation using [1(2,3)4]pentamantane shown in Fig.1. For all the simulations, we used an NPT-ensemble. Specifically, we set the initial temperature T_i = 0.12E_b ≈ 5000 K, final temperature T_f = 0.024E_b ≈ 1000K, and pressure P = 0.012E_b ≈ 10GPa. We set the C-C bond formation energy E_b = 3.6eV. In each set of the simulation for different diamondoid molecules, we kept the total number of carbon atoms approximately the same (see Table1).

For each set of parameters, we performed ten independent simulations until the newly formed C-C bond numbers are converged. Figure2a, b visualize an representative structural update for increasing bond numbers before and after one Monte Carlo step. Figure2c shows the C-C bond number as a function of theMonte Carlo step for adamantane of total 5000 carbon atoms on ten independent Monte Carlo simulations. The C-C bond number for each Monte Carlo run follows a trend, in which the bond number increases quickly at the initial stage during heating, after which it slows down and saturates, suggesting the convergence of structural optimizations. Figures1 and 3 show a few representative examples of the final diamond structures formed from diamondoid molecules. The C-C bond number convergence for each diamondoid molecule and different total numbers of carbon atoms are shown in Fig.S1.

Once the structures have converged, we characterized the final structure’s crystalline quality by first counting the number of vacancies (see Fig.4a). Here, vacancies are defined as empty carbon sites inside the final bulk diamond-like structure. We find that for #C = 10,000, the number of vacancies decreases as we move from lower diamondoids to higher diamondoids; and [1(2,3)4]pentamantane forms fewer vacancies compared with [1212]pentamantane.

To further assess the diamondoid molecules’ ability to grow into larger diamond crystals, we calculated the ratio of compatible surface sites. We defined the compatibility as the ability for a diamondoid molecule to attach to a surface site in a manner that is compatible with the diamond crystal structure. The compatible ratio is obtained by dividing the number of compatible sites to the total number of surface sites. In Fig.4b, we see that such a ratio is higher for higher diamondoids, and [1(2,3)4]pentamantane is better than [1212]pentamantane. A higher ratio means a higher probability to add diamondoids onto the surface sites, thus a faster growth rate. We note that the ratio decreases for all diamondoids as we increase the total number of carbon atoms, and when the ratio reaches zero, the diamond structure cannot grow any further. These results suggest that higher diamondoids can grow faster and form larger diamonds with fewer defects in the bulk.

Our Monte Carlo simulations show that for a fixed number of total carbon atoms, higher diamondoid molecules form a diamond crystal with fewer vacancies inside as-formed structures, and more compatible sites on the surface to potentially form a larger diamond crystal, compared to lower diamondoid molecules. Both findings suggest that higher diamondoids have a higher efficacy in forming a diamond crystal than lower diamondoids.

Our theoretical trendsalign consistently withthe experimental results of diamond crystal formation from lower diamondoids using the HPHT method², where diamond could be more easily formed from triamantane than from adamantane and diamantane, with alower critical temperature at a givenpressure. Such consistency may be attributed to the factthat the internal C-C bonds of lower diamondoids have a lowerlikelihood ofdissociating at an applied high pressure and high temperature. In diamond synthesis using the HPHT method with adamantane, ab initio molecular dynamics simulations have indicated that part ofthe C-C skeleton canremain intact during the transformation to diamond². However, whether the C-C skeleton will predominantly remain intact for larger diamondoids during diamond formation using the HPHT method has not been investigated. To further address this issue, theories utilizing large-scale molecular dynamics simulations¹⁰ may be required, where the internal carbon-carbon bond breaking is properly integrated. The different results of diamond growth using diamondoids under HPHT and CVD⁷ conditionsfurther emphasize that crystallization during these two methods intrinsically follows very different microscopic mechanisms; thus finding an optimal crystal growth method can be critical to materials synthesis and advanced materials design. Fundamental understanding of these dynamical processes associated with the CVD and HPHT diamond synthesis from diamondoids can further guide the bottom-up synthesis of diamond with desirable defect centers from doped molecular precursors^{11,12,13,14,15}. The numerical approach we present here can potentailly beapplied to a broad variety of systems for materials synthesis, such as self-assembly, HPHT and CVD processes, and more.

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 epitaxy¹⁶. 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 crystals¹⁷. 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 crystals¹⁸ and evolutionary algorithms for a broad range of materials systems^19,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 theC-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 molecules². 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←(x_max − x_min) × (y_max − y_min) × (z_max − z_min), where, for example, x_max 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 M₀ such that \({M}_{0}\cap {{{{\mathcal{S}}}}}_{0}\ne \varnothing\) and remove it.

n₀← number of bonds broken when removing M₀

\({{{{\mathcal{E}}}}}_{nn}\leftarrow\) empty nearest neighbor sites of the cluster after M₀ 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 M₀.

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 M₀ is chosen to be put at C₁ and a random orientation of M₀ around C₁ is chosen. The possible orientations are restricted to those such that the C-C bonds of M₀ align with those of the diamond lattice of the cluster).

repeat

Randomly pick an atom C₀ ∈ M₀.

Try to put C₀ at the position of C₁ and form bonds between C₀ and the existing atoms in the cluster

if Bonds can form s.t. M₀ 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

n₁← number of bonds formed after putting M₀ to the new position

ΔN_b←n₁ − n₀

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

M₀ jump to the new position

end if

end if

if M₀ does not jump to a new position then

put M₀ 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 = -E_bΔN_b + PΔV before and after the movement of this diamondoid molecule, where ΔN_b is the number of net carbon-carbon bonds formed between diamondoid molecules (if fewer C-C bonds formed than the C-C bonds breaking, ΔN_b would be negative), E_b 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/k_BT. 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 T_i, compared to the C-C bond dissociation energy. The temperature was then gradually reduced to a final low temperature T_f, at which we performed enough Monte Carlo simulations until the structure was converged. Specifically, 2 million Monte Carlo steps are performed at T_i, then temperature is linearly lowered to T_f in 2 million Monte Carlo steps. Finally, 10 million Monte Carlo steps are performed at T_f.

The source code for our Monte Carlo simulations for diamondoids is accessible at github repository: https://github.com/taaatang/Diamondoid.git.

Authors and Affiliations

1. Department of Applied Physics, Stanford University, 348 Via Pueblo, Stanford, 94305, CA, USA

   Ta Tang

2. Department of Earth and Planetary Sciences, Stanford University, 367 Panama Mall, Stanford, 94305, CA, USA

   Sulgiye Park

3. Department of Materials Science and Engineering, Stanford University, 496 Lomita Mall, Stanford, 94305, CA, USA

   Thomas Peter Devereaux

4. Stanford Institute for Materials and Energy Sciences (SIMES), SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, 94025, CA, USA

   Thomas Peter Devereaux&Yu Lin

5. Department of Physics and Quantum Theory Project, University of Florida, 2001 Museum Road, Gainesville, 32611, FL, USA

   Chunjing Jia

Contributions

C. Jia, Y. Lin and T. Tang designed the project. T. Tang developed the codes and ran the simulations. C. Jia, Y. Lin, T. Tang and S. Park made the figures. T. Tang, S. Park, T. P. Devereaux, C. Jia, and Y. Lin contributed to the writing of the manuscript.

Corresponding authors

Correspondence to Yu Lin or Chunjing Jia.

Competing interests

The authors declare no competing interests.

Peer review information

Communications Chemistry thanks Qiang Zhu and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Apeer review file is available.

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions