Simulations of nanocrystals under pressure: Combining electronic enthalpy and linear-scaling density-functional theory

In recent years, the study of nanomaterials under pressure has acquired increased scientific and technological importance.¹ In part due to their large ratio of surface to volume atoms, nanocrystalsdisplay a host of properties that differ from those of their bulk counterparts.² New dimensions are added to phase diagrams when the sizes, surface reconstructions, and terminations of the nanocrystals are taken into account.³ This has generated particular interest in nanocrystals displaying quantum confinement with diverse applications ranging from biomarkers to quantum transistors.^4–6 There is great technological potential in the possibility of using pressure to tune the physical properties of semiconducting nanocrystals, that can depend sensitively on their structures.⁴ Attaining such control at the nanoscale holds the promise of novel technological applications such as tunable photovoltaic devices,shock-absorbers,⁷ and nanoscale stress sensors.^8,9 The additional surface effects also open the door to transformation pathways that are not available to the bulk material, potentially allowing a system to become trapped in metastable states with novel properties.^10–12 Moreover, sufficiently small nanocrystals can be synthesized with few or no defects and are thus ideal models to study the kinetics of solid-solid phase transitions.^13,14 Recently, progress has been made in directly observing structural transformations in nanocrystals¹⁵ and bulk single crystals.¹⁶ Direct monitoring of transformation pathways in nanosystems is,however, still challenging with the resolution of existing experimental probes and understanding can thus greatly benefit from the insights that computer simulations provide. The nanocrystals of interest here are of an intermediate size: larger than molecules but too small to be treated satisfactorily with macroscopic concepts such as strain and stress fields in continuum models. An atomistic treatment is crucial to capture the details of the structural changes, including the shape and surface effects.

While empirical potentials are good for modelling a variety of materials, the complex bonding rearrangements associated with structural transformations of materials such as covalent semiconductors mean that ab initio methods such as density-functional theory(DFT) are essential to capture the details of the structure and dynamics with accuracy. However, the large length- and time-scales associated with the structural transformations of experimentally relevant systems pose a significant computational challenge. The

$O(N3)$ scaling of the computational effort in traditional methods such as the plane-wave pseudopotential (PWPP) formulation of DFT limits the number of atoms N that can be simulated to a few hundred and thereby seriously constrains the attainable sizes of nanocrystals. This can be addressed by working with a linear-scaling DFT code such as ONETEP,¹⁷ for which the favorable balance of cost and accuracy allows the investigation of nanocrystals with many thousands of atoms.^18,19 Even then, the challenge persists of modelling the pressure transmission between solvent molecules and nanocrystals—in analogy to experiments where nanocrystals are dissolved and placed under pressure in a diamond anvil cell. The many degrees of freedom comprising realistic solvents and the many solvent-nanocrystal collisions that need to be averaged over to sample the appropriate thermodynamic ensemble exclude a full ab initio treatment. One approach to tackle this challenge is to retain an explicit description of the solvent by embedding an ab initio simulation of the nanocrystals within a cheaper classical description of the solvent.^20–26 However, sampling rare events such as structural transformations happening over long time-scales, whilst retaining a sufficiently short time step to describe the solvent-nanocrystal collisions, generally requires unfeasibly large numbers of molecular dynamics(MD) steps to be performed. Transformations can be obtained within shorter simulation times by over-pressurizing the systems but comparability with experiment is hindered in the process. Approaches exist to surmount this issue by accelerating the free energy landscape exploration and have been applied to the pressure-induced structural transformations of nanocrystals^27,28and bulk crystals.^29–32 In practice,however, these remain computationally demanding.

Alternatively, constant pressure simulations of finite systems can be performed, in both MD and quasistatic geometry optimization, by directly optimizing the enthalpy once a suitable definition for the finite volume has been made. This can be done in a variety of ways in terms of atomic or electronic coordinates leading to an implicit description of the solvent. Some examples of total volume definitions have been suggested in terms of: a sum of atomic volumes,³³ a function of the average inter-particle distance,³⁴ the inertia tensor eigenvalues,^35,36 and the smallest convex polyhedron to circumscribe all surface atoms.³⁷These different approaches have been compared elsewhere and were shown to qualitatively reproduce results obtained with explicit solvents.³⁸ By working with quasistatic geometry optimizations at zero temperature, one removes the need for equilibration with barostats and thermostats thereby giving a comparatively inexpensive way of sampling the enthalpy landscape. Depending on system complexity, this may not give a globally optimized structure nor precise information on transition paths; however, it provides the structure and energetics of the nearest local minimum.

In the present work we use an electronic enthalpy functional H = U +P_inV_e, where U is the total Kohn-Sham internal energy of the system, P_in is the input pressure, and V_e is a volume definition based on an electronic-density isosurface.⁷ The latter allows for the description of complex geometries and the enthalpy is optimized within the linear-scaling DFT code ONETEP. We introduce an approach to calibrate the parameters defining V_e in the context of geometry optimizations and use it to simulate pressure-induced structural transformations in hydrogenated Si nanocrystals. Our results are comparable to those obtained with other methods^7,22,24and validate our approach. Si nanocrystals are of intrinsic interest due to their potential to overcome the indirect character of the lowest-energy interband transition and to be useful in optoelectronic devices.^39–41 Recently, Si nanocrystals with structures based on high-pressure bulk phases have been proposed as candidates for photovoltaic applications as they display multi-exciton generation and high quantum efficiencies.⁴² They are also found to transform under pressure between a variety of crystalline and amorphousstructures that are still the subject of theoretical and experimental investigations in both the porous and colloidal forms. Si₁₈₁H₁₁₀, the largest nanocrystal in the present work, of diameter 2.2 nm, is comparable in size to experimentally tested organically passivated colloidal nanocrystals⁴³ and demonstrates the capability of our approach.

The linear-scaling DFTcode ONETEP is based on the single particle density-matrix (DM) n(r, r′) formulation of the Kohn-Sham equations in terms of a set of local orbitals{ϕ_α(r)}, referred to as non-orthogonal generalized Wannier functions(NGWFs). These are spatially localized within spheres of radii {R_α}centered on the atomic coordinates as

where K^αβ is called the density-kernel. The electronic density ρ(r) is related to the DM by ρ(r) =2n(r, r), where the factor of two accounts for the spin degeneracy. The NGWFs are themselves expanded in terms of a fixed underlying basis of psinc functions equivalent to a systematic plane-wave basis.⁴⁴ In the course of a calculation the total energy is minimized with respect to both K^αβ and {ϕ_α} in two nested loops, subject to the constraints of normalization and idempotency.⁴⁵Linear scaling is achieved by exploiting the property of nearsightedness that allows the DM to be truncated for systems with an energy gap.⁴⁶ The electronic structure can then be described with plane-wave accuracy in terms of a minimal basis of in situoptimized NGWFs.

It has been shown that the elastic properties of bulk Si in the diamond phase calculated with ONETEP and the PWPP code CASTEP⁴⁷ give equivalent results⁴⁸ when using the same norm-conserving Si pseudopotential,⁴⁹ local density approximation exchange-correlation functional,^50,51 and plane-wave cutoff E_c. Beyond the fact that it is well-described by DFT, Si was used as a test system here due to the plethora of experimental data and computational studies with different pressure methods available for comparison.

For the calibration we require a reference DFT bulk modulus at zero pressure B₀ and its pressure derivative

$B0′=∂B/∂P|P=0$⁠. This was calculated with CASTEP using a 2-atom primitive simulation cell, a grid of 8 × 8 × 8 k-points and E_c = 800 eV. By fitting the universal Vinet equation of state^52,53 we obtained values of B₀ = 96.85 GPa and

$B_{0}^{\prime}=4.08$

$B0′=4.08$⁠.

The local orbital approach has the advantage that the {ϕ_α} are strictly zero on all grid points outside the localization radii⁵⁴ and vacuum comes at a negligible computational overhead. This is particularly advantageous for finite systems as interaction with periodic images can easily be eliminated.⁵⁵ Pulay corrections to the Hellmann-Feynman forces are required to achieve accurate ionic forces and optimized structures.^56,57 The use of in situ optimized orbitals reduces the egg-box effect⁵⁸ observed in fixed orbital approaches.

The nanocrystals were quasistatically relaxed at different pressures using the quasi-Newton Broyden-Fletcher-Goldfarb-Shanno algorithm for geometry optimization.⁵⁹ The parameters which control the accuracy of the geometry optimization must be carefully chosen for the calculations to be converged and the structures correctly relaxed. Unless specified otherwise we used E_c = 800 eV, a universal NGWF radius R_ϕ = 8 a₀, an atomic displacement tolerance of 10⁻² a₀, an energy tolerance per atom of 2 ×10⁻⁵ Ha, and a force tolerance of

$10−3 Ha a0−1$⁠.

The electronic enthalpymethod to simulate finite systems under external pressure proposed by Cococcioni et al.⁷ introduces a thermodynamic functional H = U +P_inV_e which can be minimized self-consistently within DFTalgorithms. V_e is defined as the interior of an electronic-density isosurface at a chosen cutoff density α. Introducing the Heaviside step function in terms of density values θ(ρ), V_e is calculated as

For computational purposes, the step function can be approximated by the complementary error function as

The parameter σ adjusts the sharpness of the step function and plays an important role for numerical reasons. The resulting potential contribution is

Since the potential does not explicitly depend on the nuclear positions, the compression is implicitly transmitted to the nuclei by virtue of the forces obtained from the Hellmann-Feynman theorem.⁶⁰ This can be related to the effect of Φ_V, which for a decaying density profile as in Fig. 1, is a distorted Gaussian in real space and favors the compression of the electronic density for positive pressures. The shape of Φ_V is determined by the pair of input parameters α and σ,and approximates the solvent-nanocrystal interaction. α defines the excluded volume of the solvent molecules and σcontrols the range and intensity of interaction in a manner reminiscent of soft-sphere potentials. While the method describes the solvent implicitly, providing a homogeneous and time-averaged description, the emphasis is laid on the role played by electrons as pressure mediators with an account of the shape of the nanocrystal as the pressure is applied normal to the isosurfaces. This results in a natural description that allows the seamless modelling of the excluded volume of intricate nanocrystal geometries. It also removes the need for equilibration with barostats and focuses the computational effort on the electronic structure of the nanocrystal.

Figure 2 shows the isosurface bounding V_e for a range of values of α. For larger values of α, the isosurface describes voids inside the nanocrystal, revealed by changes in slope in the plot, and results in pressure being induced internally which compensates the applied pressure. In order to describe a realistic solvent, α has to be chosen sufficiently small to apply a homogeneous compression without describing the rugosities of the nanocrystal too closely. When going to very small values of α, unphysically large excluded volumes are obtained. Figure 1 also shows that for a given choice of

$α=3×10−4a0−3$⁠,

$\sigma =5\times 10^{-4}a_0^{-3}$

$σ=5×10−4a0−3$ leads to a Φ_V(r) that clearly fails to vanish at large radii (Fig. 1(c)). From Eq. (4) it is evident that far from the nanocrystal where ρ → 0, the potentialΦ_V ∼ (P_in/σ)exp ( −α²/2σ²), and therefore a sufficiently small value of σ/α, must be chosen for the excluded volume to be well-defined. However, a sufficiently large value of σ must be chosen for the potential Φ_V(r) to be accurately integrated on the underlying real-space grid which has a spacing of Δ = 0.25a₀ in the present work. The above considerations give us a range of sensible values for α and σ, but within this range physical properties still depend on the chosen values. An approach is still needed to better resolve these depending on the system.

In principle, if the parameters α and σ defining a physical V_e were chosen correctly, an effective pressure P_eff equal to the chosen input pressure P_in would be felt within the nanocrystal.V_e could be determined for different solvent-nanocrystal interfaces by comparison with simulations using explicit solvents, e.g., with MD. This would however not be practical in a fully ab initio way as explained in Sec. I. An empirical parameterization would also be difficult considering the limited resolutions of experimental methods. Alternatively, α and σ can be calibrated by comparing P_eff and P_in if a satisfactory definition of P_eff is available. This can be done by exploiting the virial theorem in MD simulations,³⁶ but not in geometry optimization calculations. For these, a promising approach is to exploit the experimental⁴³ and computational⁶¹ result that bonds in the bulk-like core of sufficiently large alkane-terminated diamond phase Si nanocrystals display elastic properties similar to the bulk for a range of sizes. P_eff can then be estimated for a range of systems and pressures from the compression of core bonds after quasistatic geometry optimization.

Here we use the Vinet equation, with the bulk values of B₀ and

$B0′$ obtained as discussed in Sec. II, expressed in terms of the compressed and equilibrium (0 GPa) bulk-like conventional lattice parameters a and a_eq:

Defining a for the nanocrystal in terms of averaged bond lengths for core Si atoms (chosen as those atoms that are bonded exclusively to other Si atoms), an effective pressure P_eff experienced by the nanocrystal can be estimated from the volume derivative of the Vinet equation of state:

The Vinet equation holds in the absence of phase transitions and was found to give similar, albeit better fitted, results than the Birch-Murnaghan equation.⁶² P_eff can then be compared to the input pressure P_in that generated the compression, thus allowing the calibration of αand σ.

A hydrogenated tetrahedral Si₇₁H₆₀nanocrystal in the diamond phase was used for the calibration as it displays a sizable core which behaves elastically like the bulk and justifies our use of bulk values for B₀ and

$B0′$ in the calibration. The average Si–Si bond length is found to be contracted compared to the bulk; the contraction is reduced and tends towards the bulk value as the size of the nanocrystal is increased. Looking at individual Si–Si bonds, it is found that the outer shell is contracted, which has been interpreted as due to surface stress,^63,64 while inner shells substantially agree with bulk values. Similar results have been found for Si₂₉H₃₆ and Si₃₅H₃₆ using DFT with norm-conserving pseudopotentials and a local density approximation exchange-correlation functional.⁶¹ In classical linear elasticity, inhomogeneities, whether in vacuo or embedded in a material, have size-independent elastic fields.⁶⁵ This is the result of neglecting surface energies which can be justified when the ratio of surface to volume atoms is small. At the nanoscale, however, this ratio becomes important and the surface induces size-dependent elastic fields that are long-range.⁶⁶ One would expect the surfaces to induce a size-dependent strain-field and to distort the core atoms for the nanocrystal sizes investigated in this work. However, the stiffness of Si nanocrystals and the absence of reconstruction of the hydrogen-passivated surfacesresult in distortions that are smaller than the displacement tolerance. This limits the effect of the surfaces for the sizes considered and simplifies the mapping between effective pressure and compression. While the bond distributions change with the selection of the core atoms, the positions of peaks of the distribution are found to be insensitive, within the displacement tolerance, to that choice when excluding the outer shell atoms.

Figure 3 shows the results of our calibration for a range of parameter choices. A procedure where B₀ and

$B0′$ were fitted separately from Eq. (5) for each α (diamonds) is seen to produce poor agreement between P_eff and P_in. By contrast, it was found that using fixed values of B₀ and

$B_0^{\prime }$

$B0′$⁠, either from DFT bulk values from CASTEP (crosses) or experimental values (squares),produced very similar results, and the expected linear relationship was observed (Fig. 3(a)). This suggests that the assumptions entering our calibration approach and the volume definition are valid. Figure 3(b)shows that for a fixed value of α, the P_eff curves converge as a function of σ, towards the

$\sigma =5\times 10^{-5}a_0^{-3}$

$σ=5×10−5a0−3$ line. This highlights the importance of tuning σ for an accurate definition of V_e as discussed in Sec. II. Finally, given an appropriate choice of σ, Fig. 3(c) shows that there is a dependence of the compression on the chosen α. This can be understood from the volume definition of Eq. (2):changing α corresponds to using a different model of solvent-nanocrystal interaction, by altering the electronic density up to which solvent molecules approach the nanocrystal. The range of α values 3.0 × 10⁻⁴ to

$1.0\times 10^{-3}a_0^{-3}$

$1.0×10−3a0−3$ is found to give agreement between P_effand P_in within 15%.

While the parameters were calibrated on the Si₇₁H₆₀nanocrystal and model a representative solvent-nanocrystal interface, the parameters are expected to be transferable to siliconnanocrystals of different sizes and shapes, particularly if they have similar surface facets, surface reconstructions, or similar ligands. We illustrate this by comparing the effective pressure P_eff obtained at a representative value of P_in = 10 GPa, with

$α=3.0×10−4a0−3$ and

$\sigma =5.0\times 10^{-5}a_0^{-3}$

$σ=5.0×10−5a0−3$⁠, for three different nanocrystal sizes: Si₃₅H₃₆,Si₇₁H₆₀, and Si₁₈₁H₁₁₀. We obtain P_eff = 11.9, 11.5, and 11.5 GPa, respectively, indicating that the bulk-like cores of Si₇₁H₆₀ and Si₁₈₁H₁₁₀ have identical responses under pressure while Si₃₅H₃₆ displays a small discrepancy, which can be understood as due to the small nanocrystal size and consequent small region of bulk-like core. The transferability is expected to hold beyond hydrogenated silicon because all semiconductornanocrystalsurfaces, with or without organic ligands, will display a fairly similar range of values of valence electronic density inside and a similar exponential decay rate outside the surface.⁸² The solvent-nanocrystal interface resulting from a given choice of parameters can be considered appropriate so long as the lengthscale of the isosurface variations is smaller than the size of the solvent molecules given by their van der Waals surface.

We now turn our attention to structural transformations in the Si nanocrystals Si₃₅H₃₆ and Si₇₁H₆₀, which have been studied by other methods.^20–24 We then demonstrate capability for larger system sizes by studying Si₁₈₁H₁₁₀. Bulk Si is of great technological importance and has been widely used as semiconductor in both its crystalline and amorphous forms. Its phase diagram has been extensively studied and, under pressure, bulk Si has been observed to transform between 12 different structures. It transforms from the cubic diamond to the β-Sn phase at 11.7 GPa, followed by the Imma phase, primitive hexagonal (ph), orthorhombic, hexagonal closed packed (hcp), and face centered cubic (fcc) phases at respectively 13.2, 15.4, 38, 42, and 79 GPa.^67–69 Upon pressure release the ph and β-Sn phases are observed, but the diamond phase is not recovered upon full release of the pressure. Instead, a host of crystalline and amorphous metastable structures are observed, the most common of which are the BC8 and R8 phases that correspond to distorted tetrahedral structures. Unlike the bulk, the small size and the stabilising effect of the surfaces of Si nanocrystals allows for metastable structures and transformation mechanisms that are still the subject of investigation. Size-dependence of structural, optical, and electronic properties has been reported for a range of nanocrystal sizes in both colloidal^3,43 and porous forms.^70,71

Recent X-ray diffraction (XRD) experiments on colloidal plasma-synthesized Si nanocrystals, where the surfaces are initially H terminated and later functionalized with dodecane, investigate nanocrystals of diameters 3.2, 3.8, and 4.5 nm, under pressures in the range 0–73 GPa at room temperature.⁴³ A transformation between the diamond and what is interpreted to be the ph structure is reported in the range 17–22 GPa although the small size of the sample results in significant broadening of the spectra and makes it difficult to identify the structure. Another structural transformation occurs in the range 40–44 GPa and matches an hcp structure. The ph structure is recovered upon decompression down to as low as 18.4 GPa followed by a stable amorphousstructure upon complete decompression.

XRD and Raman spectroscopy experiments⁷⁰performed on porous Si with average crystallite diameters ∼5 nm (with distributions at 3 nm and 7 nm) report a transformation from diamond to a high-density amorphous (HDA) phase at 14 GPa which, upon pressure release, transforms to a low-density amorphous(LDA). More recent work⁷¹ on porous Si with crystallites of ∼4 nm average diameter observe a transformation from diamond to ph phase at 20 GPa with no amorphisation up to 39 GPa. A HDA phase is recovered upon decompression around 15 GPa followed by an LDA phase at 4.5 GPa. Under a further pressurization cycle, an amorphous to crystalline transformation is observed between LDA and ph at 18 GPa. Such reversible transformations between LDA, HDA, and ph phases have also been observed for amorphous bulk Si (a-Si).^72–74 

Theoretical investigations have attempted to characterize these amorphous phases in bulk Si^73–78 and hydrogenated Si nanocrystals.^20–24 The LDA phase has been described as a disordered tetrahedrally coordinated network, the HDA as a deformed tetrahedral network⁷⁵ with the presence of interstitials increasing the coordination to 5–6, and finally a very-high-density amorphous phase (VHDA) has been postulated⁷⁸ with coordinations 8–9 as found also in ice.⁷⁹ This classification of the amorphousstructures is adopted in this paper.

Structural properties are generally analysed in terms of bond-length distribution, coordination number, and bond-angle distribution. Moreover, we employ ring statistics as a way of tracking the evolution of the covalent Si networks. Every Si atom can be treated as a node and bonds as links connecting the nodes. We define an n-membered ring as a closed path connecting n atoms according to Guttman's definition, focusing on the total number of rings R_T and the proportion of nodes belonging to at least one n-membered ring P_n.^80,81 Two atoms are considered to be bonded when separated by a distance smaller than the first minimum of the radial distribution function of the bulk R_c = 5.33a₀. All nanocrystals were initially relaxed with geometry optimization at 0 GPa, and then pressure was applied in steps of 5 GPa or less, relaxing the geometry at every pressure to find the minimal enthalpy configuration. Si₃₅H₃₆ and Si₇₁H₆₀ were investigated in a pressure range 0–50 GPa while only 0–25 GPa was considered for Si₁₈₁H₁₁₀ as the system becomes metallic beyond 25 GPa and occupancy smearing would be needed.

A density cutoff value of

$α=3.0×10−4a0−3$ was chosen for the rest of the calculations because it was shown to result in a good calibration (Fig. 3) and as it enables direct comparison with the simulations of Cococcioni et al.⁷ 

Figure 4 shows the structures of Si₇₁H₆₀as it is compressed in the range 0–50 GPa and depressurized to 5 GPa, while Fig. 5 shows a range of descriptors of these transitions. The bond-angle distribution in Fig. 5(a) initially shows a single peak at 109.5°, typical of the tetrahedral diamondstructure. The peak broadens when the pressure applied increases from 0 GPa to 20 GPa, which reflects a lower degree of symmetry in the structure. No change in the coordination number of Si atoms is observed at this stage and a tetrahedral coordination is retained for the innermost Si atoms. When the pressure applied is increased to ∼35 GPa, a structural change is observed as the nanocrystalamorphizes, as evidenced by the broad bond-angle distribution. The transformation is mediated by the breaking and making of bonds which are accompanied by the appearance of interstitial atoms and an increase in coordination numbers (Fig. 5(b)). To better resolve the transformation, calculations for Si₇₁H₆₀ were repeated in steps of 2 GPa in the interval 20–30 GPa. Between 21 and 30 GPa, the average coordination of Si atoms (Fig. 5(c)) increases from 4 at 0–20 GPa, to 5.1 at 27 GPa which is consistent with a transformation from diamond to HDA. However, even at 30 GPa the nearest neighbour peak remains,suggesting that some local order is retained and the first coordination shell is preserved. As the pressure is increased further, the average coordination reaches 8.3 at 50 GPa which matches that of a VHDA structure. Upon pressure release, the average coordination plummets to 6.7 at 20 GPa and 4.3 at 5 GPa corresponding to a disordered tetrahedral network consistent with an LDAstructure. Similarly, for Si₃₅H₃₆ one obtains a coordination of 4 at 0 GPa, 5.5 at 30 GPa, 7.3 at 40 GPa, 7.7 at 50 GPa, 4.5 upon decompression to 5 GPa, and for Si₁₈₁H₁₁₀ 4 at 0 GPa and 5.5 at 25 GPa. Figure 5(d) shows the electronic volume per atom for Si₇₁H₆₀. This reveals discontinuous changes in the intervals 20–21 GPa and again at 27–30 GPa, which suggests that the structural transformations from diamond to HDA and HDA to VHDA are first order.

Our results are consistent with previously reported Car-Parrinello MD simulations on Si₃₅H₃₆ and Si₇₁H₆₀ at 600 K using a classical explicit soft-sphere solvent²² and on Si₃₅H₃₆ at 300 K using the electronic enthalpymethod.⁷ In the former simulations, transformations from the diamond to an amorphousstructure with average coordination of core Si atoms of 7.3 were reported at about 30 GPa for Si₇₁H₆₀and 35 GPa for Si₃₅H₃₆. An amorphousstructure with average coordination 4.3 is recovered upon decompression to 5 GPa. In the latter simulations,Si₃₅H₃₆ is found to amorphize around 26–28 GPa with average coordination of Si atoms reaching ∼6.5 and to remain amorphous upon pressure release to 0 GPa with an average coordination of ∼4. No sensitivity of the results was reported when changing α. This can be understood as resulting from the small range of values chosen, combined with thermal noise concealing the dependence on parameters seen in the present work. Differences in transformation pressure for the same system between MD simulations using an explicit solvent²² and the electronic enthalpy method⁷ are due to the different way that pressure is applied.³⁶ We do not expect direct agreement in the transformation pressure with MD nor experimental results considering that we have used a quasistatic approach: in the absence of thermal fluctuations, one needs to overpressurize the system to overcome the large activation energies associated with the energetic cost of making and breaking bonds. The absence of thermal noise in our simulations does, however, facilitate the detailed monitoring of the amorphization and provides complementary information to that obtained with MD.

The ring statistics shown in Figs. 5(e) and 5(f)indicate the presence of four distinct regions. In the interval 0–20 GPa, only 6-membered rings are present. The population of 6-membered rings decays in favor of 3-, 4- and 5-membered rings in the interval 21–30 GPa. Between 30 and 50 GPa, the 3-membered ring population grows further at the expense of 4- and 5-membered rings and dominates at 50 GPa. Upon pressure release, the 4-, 5-, 6-,and 7-membered ring population recovers while the 3-membered ring population drops sharply as the nanocrystal dilates. 6-membered rings are a signature of the corner-sharing tetrahedra in the diamond cubic structure, while 3-membered rings arise through the formation of the equilateral triangles that cause the peak at 60° in the bond-angle distribution (Fig. 5(a)). The presence of 3-, 4-, and 5-membered rings indicates amorphization and the larger 7-membered ring the formation of voids (Fig. 6). Our results are consistent with the existence of three amorphousstructures visited during the pressure-induced structural transformation: HDA corresponding to average coordination numbers ∼5–6,VHDA with coordinations ∼8–9, and LDA with coordination ∼4 obtained upon pressure release.

A reversible amorphization diamond→HDA→diamond is obtained when performing a pressure cycle between 0 and 30 GPa (Fig. 7). By contrast, when increasing the pressure all the way to 50 GPa, irreversible bonding events accompany the transformation, which proceeds as diamond→VHDA→LDA. The final LDAstructure is found to be higher in energy compared to the original diamondstructure and thus corresponds to a local minimum in energy. The nanocrystals display hysteresis and comparing the upstroke and downstroke 20 GPa structures, it is clear that they are respectively crystalline and amorphous. This behavior demonstrates the possibility of trapping nanocrystals in unconventional bonding geometries when performing a pressure cycle,yielding electronic properties different from the bulk. Of great promise is the possibility of designing materials with desired properties using pressure as a way to tune these.

In particular, it is found that the HOMO–LUMO gap changes dramatically under pressure and the way this happens is strongly size-dependent. While it is well known that the local density approximation underestimates the size of the gap for Si, it does give qualitative information and significant trends. From Fig. 8 we observe that beyond a certain pressure,the gap drops sharply to a lower value. All clusters are semiconductors at low pressure with increased gaps for the smaller nanocrystalscompared to the bulk. This can be understood as due to quantum confinement which is significant for the nanocrystals under investigation. The Si₁₈₁H₁₁₀nanocrystal gap sharply drops to 0.04 eV at ∼25 GPa showing that it becomes metallic, whereas the smaller clusters retain a sizeable gap even at higher pressures. In the pressure range 0–20 GPa, the gap of Si₃₅H₃₆ increases with a slope that reduces as the pressure is ramped up. Meanwhile, for Si₇₁H₆₀, the gap increases at first and decreases in the range 10–20 GPa. For Si₁₈₁H₁₁₀, the gap decreases linearly in the range 5–20 GPa. This size-dependent behavior can be interpreted as a competition between quantum confinement and the negative pressure coefficient of diamond Si. The former tends to increase the gap as the nanocrystals are compressed, while the latter tends to decrease it due to the dominance of the indirect transition (corresponding to X_conduction − Γ_valence in the bulk). Si₁₈₁H₁₁₀, of diameter 2.2 nm, has a change of gap with pressure of −10.7 meV GPa⁻¹ which is of the same order as the experimental results for 2.6 nm diameter nanocrystals of −17.2 meV GPa⁻¹.⁴³ As the nanocrystal size is increased, the quantum confinement effect is expected to vanish and a linear decrease of the gap remain with a slope tending to the DFT value for bulk diamond Si of−15.4 meV GPa⁻¹. The above results highlight the capability of our approach to simulate sizes of experimental relevance with DFT accuracy.

We have implemented an electronic enthalpy method in a linear-scaling DFT code (ONETEP) to simulate nanocrystals under pressure. An approach to calibrate the parameters defining the electronic volume in the context of geometry optimizations was proposed, demonstrating how the pressure felt inside the nanocrystal can be successfully matched to the input pressure in the electronic enthalpyfunctional. We have applied this method to the structural transformations of hydrogenated Si nanocrystals of different sizes and obtained results in good agreement with simulations using explicit solvents. Our quasistatic investigation has allowed for the detailed study of polyamorphic transformations and provided information that would be difficult to extract with the thermal noise of a MD simulation. Size-dependent structural transformations were obtained between the diamondstructure and the amorphousLDA, HDA, and VHDA structures. The behavior of the intermediate structuresupon pressure release was investigated and depressurizing from HDA and VHDA structures was found to give respectively diamond and LDAstructures. These have distinct electronic properties and the changes in HOMO–LUMO gaps with pressure were analyzed for different nanocrystal sizes and display qualitative agreement with experiment of similar diameters. The present work highlights the capability of our approach; barring further progress in the synthesis and probing of smaller nanocrystal sizes,techniques such as linear-scaling DFT become essential to simulate sizes of experimental relevance with quantum accuracy.

The authors would like to thank Francesco Mauri (Université Pierre et Marie Curie, Paris, France)for a helpful discussion. The authors are grateful for the computing resources provided by the Imperial College High Performance Computing Service and the HECToR UK supercomputer facility which have enabled all the simulations presented here. The HECToR computer time was funded by the UK Engineering and Physical Sciences Research Council (EPSRC) Grant No. EP/F037457/1 and a RAP class 1b award. N.R.C.C. was supported through a studentship in the Centre for Doctoral Training on Theory and Simulation of Materials at Imperial College funded by EPSRC Grant No. EP/G036888/1. N.D.M.H. acknowledges the support of EPSRC Grant No. EP/G05567X/1, a Leverhulme early career fellowship, and the Winton Programme for the Physics of Sustainability. P.D.H. acknowledges the support of a Royal Society University Research Fellowship.