Simulations of Nanocrystals Under Pressure: Combining Electronic Enthalpy and Linear-Scaling Density-Functional Theory

We present an implementation in a linear-scaling density-functional theory code of an electronic enthalpy method, which has been found to be natural and efficient for the ab initio calculation of finite systems under hydrostatic pressure. Based on a definition of the system volume as that enclosed within an electronic density isosurface [Phys. Rev. Lett., 94, 145501 (2005)], it supports both geometry optimizations and molecular dynamics simulations. We introduce an approach for calibrating the parameters defining the volume in the context of geometry optimizations and discuss their significance. Results in good agreement with simulations using explicit solvents are obtained, validating our approach. Size-dependent pressure-induced structural transformations and variations in the energy gap of hydrogenated silicon nanocrystals are investigated, including one comparable in size to recent experiments. A detailed analysis of the polyamorphic transformations reveals three types of amorphous structures and their persistence on depressurization is assessed.


I. INTRODUCTION
In recent years, the study of nanomaterials under pressure has acquired increased scientific and technological importance. 1 In part due to their large ratio of surface to volume atoms, nanocrystals display a host of properties that differ from those of their bulk counterparts. 2New dimensions are added to phase diagrams when the sizes, surface reconstructions and terminations of the nanocrystals are taken into account. 35][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. 4Attaining such control at the nanoscale holds the promise of novel technological applications such as tunable photovoltaic devices, shock-absorbers, 7 and nanoscale stress sensors 8,9 .1][12] Moreover, sufficiently small nanocrystals can be synthesized with few or no defects and are thus ideal models to study the kinetics of solidsolid phase transitions. 13,14Recently, progress has been made in directly observing structural transformations in nanocrystals 15 and bulk single crystals 16 .Direct monitoring of transformation pathways in nanosystems is, however, still challenging with the resolution of existing a) Electronic mail: niccolo.corsini@imperial.ac.uk 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 timescales associated with the structural transformations of experimentally relevant systems pose a significant computational challenge.The O(N 3 ) scaling of the computational effort in traditional methods such as the planewave 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, 17 for which the favorable balance of cost and accuracy allows the investigation of nanocrystals with many thousands of atoms. 18,19Even 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 solventnanocrystal collisions that need to be averaged over to sample the appropriate thermodynamic ensemble exclude a full ab initio treatment.][22][23][24][25][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.][31][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, 33 a function of the average inter-particle distance, 34 the inertia tensor eigenvalues 35,36 and the smallest convex polyhedron to circumscribe all surface atoms. 37These different approaches have been compared elsewhere and were shown to qualitatively reproduce results obtained with explicit solvents. 38By 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 in V e , where U is the total Kohn-Sham internal energy of the system, P in the input pressure and V e a volume definition based on an electronic-density isosurface. 7The 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,24 and validate our approach.0][41] Recently, Si nanocrystals with structures based on high-pressure bulk phases have been proposed as candidates for photovoltaic appli-cations as they display multi-exciton generation and high quantum efficiencies. 42They are also found to transform under pressure between a variety of crystalline and amorphous structures that are still the subject of theoretical and experimental investigations in both the porous and collolloidal forms.Si 181 H 110 , the largest nanocrystal in the present work, of diameter 2.2 nm, is comparable in size to experimentally-tested organically passivated colloidal nanocrystals 43 and demonstrates the capability of our approach.

II. METHODOLOGY
The linear-scaling DFT code 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. 44In 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. 45Linear scaling is achieved by exploiting the property of nearsightedness that allows the DM to be truncated for systems with an energy gap. 46The electronic structure can then be described with plane-wave accuracy in terms of a minimal basis of in situ optimized NGWFs.It has been shown that the elastic properties of bulk Si in the diamond phase calculated with ONETEP and the PWPP code CASTEP 47 give equivalent results 48 when using the same normconserving Si pseudopotential 49 , local density approximation exchange-correlation functional 50,51 and planewave 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 0 and its pressure derivative B ′ 0 = ∂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 0 = 96.85GPa and B ′ 0 = 4.08.The local orbital approach has the advantage that the {φ α } are strictly zero on all grid points outside the local- ization radii 54 and vacuum comes at a negligible computational overhead.This is particularly advantageous for finite systems as interaction with periodic images can easily be eliminated. 55Pulay corrections to the Hellmann-Feynman forces are required to achieve accurate ionic forces and optimized structures. 56,57The use of in situ optimized orbitals reduces the egg-box effect 58 observed in fixed orbital approaches.
The nanocrystals were quasistatically relaxed at different pressures using the quasi-Newton BFGS algorithm for geometry optimization. 59The 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 0 , an atomic displacement tolerance of 10 −2 a 0 , an energy tolerance per atom of 2 × 10 −5 Ha and a force tolerance of 10 −3 Ha a −1 0 .The electronic enthalpy method to simulate finite systems under external pressure proposed by Cococcioni et al. 7 introduces a thermodynamic functional H = U +P in V e which can be minimized self-consistently within DFT algorithms.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 60 .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 favours 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 timeaveraged 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 −4 a −3 0 , σ = 5 × 10 −4 a −3 0 leads to a Φ V (r) that clearly fails to vanish at large radii (Fig. 1c).From Eq. 4 it is evident that far from the nanocrystal where ρ → 0, the potential Φ V ∼ (P in /σ) exp(−α 2 /2σ 2 ), and therefore a sufficiently small value of σ/α must therefore 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.25 a 0 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.

III. CALIBRATION
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 Section I.An empirical parametrization 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 36 , but not in geometry optimization calculations.For these, a promising approach is to exploit the experimental 43 and computational 61 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 0 and B ′ 0 obtained as discussed in Section II, expressed in terms of the compressed and equilibrium (0 GPa) bulk-like conventional lattice parameters a and a eq : a a eq − 1 . (5) 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. 62P 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 71 H 60 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 0 and B ′ 0 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 29 H 36 and Si 35 H 36 using DFT with norm-conserving pseudopotentials and a local density approximation exchange-correlation functional. 61n classical linear elasticity, inhomogeneities, whether in vacuo or embedded in a material, have size-independent elastic fields. 65This 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 sizedependent elastic fields that are long-range. 66One would expect the surfaces to induce a size-dependent strainfield 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 surfaces result 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 changes 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 0 and B ′ 0 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 0 and B ′ 0 , 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 σ = 5 × 10 −5 a −3 0 line.This highlights the importance of tuning σ for an accurate definition of V e as discussed in Section 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 −4 to 1.0 × 10 −3 a −3 0 is found to give agreement between P eff and P in within 15%.
While the parameters were calibrated on the Si 71 H 60 nanocrystal and model a representative solventnanocrystal interface, the parameters are expected to be transferable to silicon nanocrystals 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 −4 a −3 0 and σ = 5.0 × 10 −5 a −3 0 , for three different nanocrystal sizes: Si 35 H 36 , Si 71 H 60 and Si 181 H 110 .We obtain P eff = 11.9, 11.5 and 11.5 GPa respectively, in-dicating that the bulk-like cores of Si 71 H 60 and Si 181 H 110 have identical responses under pressure while Si 35 H 36 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 semiconductor nanocrystal surfaces, 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. 82The 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.

IV. RESULTS: PRESSURE-INDUCED TRANSFORMATIONS IN SILICON NANOCRYSTALS
1][22][23][24] We then demonstrate capability for larger system sizes by studying Si 181 H 110 .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.][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. 43A 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 amorphous structure upon complete decompression.
XRD and Raman spectroscopy experiments 70 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 71 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.1][22][23][24] The LDA phase has been described as a disordered tetrahedrally coordinated network, the HDA as a deformed network 75 with the presence of interstitials increasing the coordination to 5-6 and finally a very-high-density amorphous phase (VHDA) has been postulated 78 with coordinations 8-9 as found also in ice. 79This classification of the amorphous structures 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,81Two 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.33 a 0 .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 35 H 36 and Si 71 H 60 were investigated in a pressure range 0-50 GPa while only 0-25 GPa was considered for Si 181 H 110 as the system becomes metallic beyond 25 GPa and occupancy smearing would be needed.
A density cutoff value of α = 3.0×10 −4 a −3 0 was chosen for the rest of the calculations because shown to result in a good calibration (Fig. 3) and as enables direct comparison with the simulations of Cococcioni et al. 7 Figure 4 shows the structures of Si 71 H 60 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 diamond structure.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 nanocrystal amorphizes, 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 71 H 60 were repeated in steps of 2 GPa in the interval 20-30 GPa.Between 21-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 LDA structure.Similarly, for Si 35 H 36 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 181 H 110 4 at 0 GPa and 5.5 at 25 GPa.Fig. 5(d) shows the electronic volume per atom for Si 71 H 60 .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 35 H 36 and Si 71 H 60 at 600 K using a classical explicit soft-sphere solvent 22 and on Si 35 H 36 at 300 K using the electronic enthalpy method 7 .In the former simulations, transformations from the diamond to an amorphous structure with average coordination of core Si atoms of 7.3 were reported at about 30 GPa for Si 71 H 60 and 35 GPa for Si 35 H 36 .An amorphous structure with average coordination 4.3 is recovered upon decompression to 5 GPa.In the latter simulations, Si 35 H 36 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 22 and the electronic enthalpy method 7 are due to the different way that pressure is applied. 36We do not expect direct agreement in the transformation pressure with MD nor experimental results considering that we have used 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 Fig. 5(e) and Fig. 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 favour of 3-, 4-and 5-membered rings in the interval 21-30 GPa.Between 30-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 amorphous structures 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 LDA structure is found to be higher in energy compared to the original diamond structure 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 nanocrystals compared to the bulk.This can be understood as due to quantum confinement which is significant for the nanocrystals under investigation.The Si 181 H 110 nanocrystal gap sharply drops 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 35 H 36 increases with a slope that reduces as the pressure is ramped up.Meanwhile, for Si 71 H 60 , the gap increases at first and decreases in the range 10-20 GPa.For Si 181 H 110 , the gap decreases linearly in the range 5-20 GPa.This sizedependent behaviour 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 181 H 110 , of diameter 2.2 nm, has a change of gap with pressure of -10.7 meV GPa −1 which is of the same order as the experimental results for 2.6 nm diameter nanocrystals of -17.2 meV GPa −1 . 43As 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 −1 .The above results highlight the capability of our approach to simulate sizes of experimental relevance with DFT accuracy.

V. CONCLUSION
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 enthalpy functional.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.Sizedependent structural transformations were obtained between the diamond structure and the amorphous LDA, HDA and VHDA structures.The behaviour of the intermediate structures upon pressure release was investigated and depressurizing from HDA and VHDA structures was  HOMO-LUMO gaps of Si35H36, Si71H60 and Si181H110 under pressure and comparison with experiment 43,83 and DFT calculations of bulk Si in the diamond structure.
found to give respectively diamond and LDA structures.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.

FIG. 5 .
FIG. 5. Structural transformations of Si71H60 under pressure (the shaded region corresponds to the depressurization): (a) distribution of the bonded Si-Si-Si angles at 20, 25 and 35 GPa; (b) distribution of the coordination numbers of Si atoms at 20, 25 and 35 GPa; (c) total number of rings RT with pressure; (d) proportion of nodes belonging to at least one n-membered ring Pn; (e) average coordination number of Si atoms with pressure; (f) electronic volume per atom with pressure.

FIG. 7 .
FIG. 7. Distribution of Si-Si distances in Si71H60 under a pressure cycle at pressures 0, 20 and 30 GPa upstroke (top panel) and downstroke starting from 30 GPa (bottom panel).