Structure and Bonding in Amorphous Red Phosphorus

Abstract Amorphous red phosphorus (a‐P) is one of the remaining puzzling cases in the structural chemistry of the elements. Here, we elucidate the structure, stability, and chemical bonding in a‐P from first principles, combining machine‐learning and density‐functional theory (DFT) methods. We show that a‐P structures exist with a range of energies slightly higher than those of phosphorus nanorods, to which they are closely related, and that the stability of a‐P is linked to the degree of structural relaxation and medium‐range order. We thus complete the stability range of phosphorus allotropes [Angew. Chem. Int. Ed. 2014, 53, 11629] by now including the previously poorly understood amorphous phase, and we quantify the covalent and van der Waals interactions in all main phases of phosphorus. We also study the electronic densities of states, including those of hydrogenated a‐P. Beyond the present study, our structural models are expected to enable wider‐ranging first‐principles investigations—for example, of a‐P‐based battery materials.


Introduction
Phosphorus is one of the structurally most diverse elements, and its various allotropes continue to attract widespread research interest in chemistry. [1]White phosphorus is the thermodynamic standard state and consists of tetrahedral P4 molecules. [2]Black phosphorus shows a layered structure with buckled six-membered rings held together by van der Waals (vdW) dispersion interactions, [3] and can be exfoliated to form monolayer "phosphorene", which is beginning to be used for multiple advanced technologies. [4]Hittorf's violet [5] and Ruck's fibrous phosphorus [6] both contain cage-like fragments, with five-membered rings as the principal building unit, and characteristic "P8" and "P9" cages being found in both modifications.These fragments connect to form perpendicular (parallel) tubular structures in violet (fibrous) phosphorus, respectively.Similar cages and tubular networks are found in recently synthesised phosphorus nanorods [7] and nanowires. [8] addition to the crystalline allotropes, there is a widely-known non-crystalline form, namely, red amorphous phosphorus (a-P). [9]An emerging application for this material is in batteries, [10] based on its ability to form Li-P and Na-P phases that lead to high theoretical capacities of a-P-based anodes, as long as the conductivity and the volumetric change during cycling can be controlled. [11]In terms of structural chemistry, various models have been proposed for a-P, including two-dimensional structures with layered motifs, similar to those in black phosphorus; [12] tubular networks containing cage-like motifs, e.g., P8 and P9 fragments; [9b] or a structural model primarily composed of P3 rings and P4 tetrahedra that form extended chains. [13]e here write "P4" rather than "P4", for consistency of notation.)However, whilst the crystalline structures can be accurately characterised in advanced diffraction and imaging experiments, [6,14] a large part of the difficulty in studying a-P is in determining its structure in the first place.Early work using neutron diffraction suggested the existence of P8 or P9 motifs, inferred from a similarity to Hittorf-type fragments, [9b] whereas empirical potential structure refinement based on X-ray diffraction data also implied the presence of P4 tetrahedra. [13]In addition, Raman spectroscopy studies suggested the existence of both buckled six-membered rings [15] and cage-like motifs, [16] indicating a rather complex atomic structure of a-P, whose details may well depend on the synthesis conditions.
We have recently shown that machine-learning (ML) methods which are "trained" on quantummechanical reference data can lead to an unprecedented level of quality in describing a-P on the atomic scale. [17]Specifically, we created an a-P structural model (containing 1,984 atoms) by simulated slow cooling from a disordered melt. [18]The resultant structure primarily contains cluster fragments of five-membered rings, in line with the long-established Baudler rules [19] and with the foundational theoretical work on gas-phase clusters by Böcker and Häser. [20]The validity of the structural model was shown by comparison to the available experimental evidence: the simulated first sharp diffraction peak (FSDP) in the structure factor of our model matches previous experimental results very well, and so does its evolution in compression and decompression simulations -see Ref. [18] and references therein.
And yet, just like there are open questions about many crystalline phosphorus allotropes, [1,21] there remain fundamental chemical questions about the amorphous form, a-P.For example, how does the degree of local structural ordering determine the energetic (meta-) stability?What is the chemical-bonding nature, which might be expected to include both strong covalent and weaker dispersion interactions?What is the role of coordination defects?And where does a-P fall within the previously introduced first-principles stability range of the crystalline phosphorus allotropes? [21] the present work, we address precisely those questions by using a suite of state-of-the-art computational chemistry techniques: ML-driven molecular-dynamics (MD) simulations, firstprinciples DFT including many-body-dispersion (MBD) corrections, as well as analyses of electronic structure and orbital interactions.We introduce a series of structural models for a-P and also for its hydrogenated form (a-P:H) -large enough to represent the complex structures, yet small enough to enable full first-principles investigations.We derive new insight into the role of local structural order in the energetic stability of a-P, the nature of defects, and the balance between covalent and vdW dispersion interactions in phosphorus allotropes.Structures of a-P as obtained from ML-driven simulations and subsequent DFT structural optimisation.We generated three independent structural models containing 248 atoms each, using the same simulation protocol in each case.The structures are labelled 1 to 3, with an increasing degree of structural order.Colour-coding indicates the coordination numbers, viz.N = 2 (blue), N = 3 (pale blue), and N = 4 (pink), determined using a 2.4 Å cut-off.(b) Relative abundance (per-atom count) of cluster fragments in the a-P models (circles), as compared to the result of a similar but larger-scale 1,984-atom simulation in Ref. [18] (squares).Solid and dashed lines are guides to the eye.Sketches of selected fragments are shown, following Ref.[20].

Results and Discussion
We carried out melt-quench MD simulations to create structural models of a-P. [18]Our simulations start from a metastable liquid, rather than the fluid (consisting of P4 molecules) which is stable at low pressure -this is done on purpose, because the large structural diversity in the metastable liquid allows the simulation to rapidly explore relevant local structural fragments.
The protocol followed our earlier work, [18] but the new simulations were now carried out in smaller simulation cells, with the aim to enable subsequent DFT studies, and in three parallel runs.Having three small-scale models of a-P allows us to take advantage of the statistical nature of the process: with the same protocol but different starting configurations, these structures over millions of simulation steps evolved into rather distinct disordered networks, which we label 1 to 3 (Figure 1).Our a-P models, after DFT-based structural optimisation, contain mostly three-fold coordinated atoms, in line with standard valence rules; there are only ≈ 1% of over-coordinated (N = 4) and ≈ 2% of under-coordinated atoms (N = 2).There is a different degree of structural ordering in the three a-P models (Figure 1a): 1 shows a more disordered, random network, whereas the arrangement of fragments in 3 follows patterns resembling those in violet and fibrous phosphorus (viz.tubular networks).Similar to the previously simulated large-scale a-P model, [18] the three models created in this work contain abundant five-membered rings, and more complex cluster fragments that are made up of those (e.g., the P8 cage found in violet and fibrous P consists of four fused five-membered rings).We use the same notation as in Ref. [20] to label those fragments.
The distribution of cluster fragments (Figure 1b), together with a similarity in short-range structural features (Figure S1c-d), confirms that our new small-scale models are overall consistent with the previously validated, larger-scale model.A comparatively higher count of P8 and P9 cages in 3 implies a greater similarity with the closely related crystalline forms, violet and fibrous phosphorus.Such cage-like fragments have also been used as building blocks of predicted phosphorus allotropes with more complex local environments, via a random structure search. [22] P4 units are observed in 2 and 3, implying pure network, rather than molecular, structures.
The overall agreement between the present structures and the one from Ref. [18] underlines that purpose-tailored simulations are possible with ML potentials: enabling ultra-large-scale simulations, but also the use of smaller-scale models (e.g., 1-3 in this work) which are accessible to subsequent first-principles DFT analyses of energetics and chemical bonding.To compare the energetic stability of different phosphorus phases, including our a-P structures (1 to 3) and relevant crystalline forms, we computed the energies of fully DFT-optimised structures for those modifications (Figure 2).Corrections for vdW interactions are needed to accurately describe longer-range interactions in phosphorus modifications: [17,23] we here used the "D3" method [24] for structural optimisations, and the many-body dispersion (MBD) method [25] for subsequent single-point computations.Both methods describe the exfoliation of black phosphorus remarkably well, [26] and MBD also captures the unit-cell volume of black phosphorus in almost quantitative agreement with experiment, and the exfoliation energy as compared to higher-level quantum-chemistry methods, [17] viz.quantum Monte Carlo [27] and coupled-cluster theory. [28]e computed energies provide a full stability range of phosphorus allotropes.Among the phases studied, white phosphorus is the least energetically stable (Figure 2).Violet, fibrous, and black phosphorus are all predicted to have close energies, with differences less than 0.3 kJ mol -1 , although the structure of black phosphorus differs from that of the other two.In fact, such energetic near-degeneracy is not only seen with HSE06+MBD, but was reported using various DFT methods. [21,29] ollowing Ref. [21], we also generated three models of P nanorods (referred to as n1, n2, and n3), by removing the copper and iodine atoms from (CuI)8P12, (CuI)3P12, and (CuI)2P14, respectively, and then relaxing the structures using DFT.These models consist of more complex tubular chains and are energetically less stable than the "textbook" crystalline allotropes.Our a-P models (1 to 3), energetically, sit in between white P and the other crystalline forms, and with increased structural ordering (Figures 1a and S2), the energies of these three models decrease.Table 1 shows that the computed mass densities of different allotropes agree well with experimental data, and the computed energy ranking of the crystalline forms is consistent with previous studies [21,23] -our combined ML-and DFT-based simulation approach has now allowed us to extend this ranking to the previously poorly understood amorphous form.
Bachhuber et al. have previously studied dispersion interactions in phosphorus allotropes, using pairwise dispersion-correction methods available at the time. [32]With advanced MBD corrections, we also observe that the vdW contributions to the overall energy vary across the different phosphorus modifications (Table 1): white phosphorus shows the smallest absolute vdW contribution (lowering the total energy by about 14 kJ mol -1 ), whereas black phosphorus has the largest value (≈ 20 kJ mol -1 ).All structures studied that contain cage-like motifs have a similar level of vdW contributions, of about 16-17 kJ mol -1 in stabilisation relative to the pure, uncorrected HSE06-DFT energy.These results can be understood from the different interspaces of building fragments in molecular and network solids: in white phosphorus, the P4 tetrahedra are well-separated from each other, and the distance between building units is longer than in other modifications, resulting in smaller dispersion contributions.This result is intuitive in the sense that vdW forces approximately decay with the sixth power of the distance, but it is still notable that the "molecular" solid, P4, is less strongly vdW bonded than any of the black, violet, fibrous, or amorphous forms.
Table 1.Properties of the a-P models 1 to 3 and comparison to crystalline phases.The standard phase, white phosphorus, was set as the reference.The relative energies of other structures are given with respect to white phosphorus.
White (β-P4) ± 0 (reference) -13.7 1.89 1.98 [2] 1 -8.8 -16.4 2.25 2.14 to 2.34 [13,30] [31] We next created structural models for hydrogenated a-P (referred to as "a-P:H"), aiming to better understand atomic coordination defects (i.e., over-and under-coordinated atoms) and how they might affect the electronic properties of a-P.Comparable hydrogenated phases of amorphous silicon (known as "a-Si:H") have long been studied: the material can be generated under synthesis conditions (e.g., plasma decomposition of silane) in which atomic defects are passivated by hydrogen atoms, [33] leading to varying concentrations of hydrogen and defects in the material. [34]Experiments have shed light on the activation and functionalisation of white phosphorus [35] to form different P-containing species (e.g., organophosphorus compounds); [36] under specific synthesis conditions, P4 tetrahedra break and give rise to unsaturated phosphorus atoms, which readily form new bonds.We here considered hydrogenated models for a-P, following a similar idea as in a-Si:H, as unsaturated (N = 2) phosphorus atoms may also be expected to be passivated during synthesis (e.g., when a-P is formed by thermal decomposition of PH3). [37]gure 3a illustrates the hydrogenation process, starting from pristine a-P (models 1 to 3), and leading to three hydrogenated structures which we label 1H to 3H.We generated these structures as follows: (1) all fourfold connected atoms (≈ 1% of the total atoms in the model) were removed; (2) the valences of all the resulting twofold connected atoms were saturated by adding one additional hydrogen atom to each, forming P-H bonds perpendicular to the two P-P bonds and on the side with more open space; (3) the modified structures were further relaxed using dispersion-corrected DFT (PBE+D3; Methods section).This procedure leaves the connectivity of different cluster fragments intact, whilst ensuring that all P atoms are three-fold connected in the resulting hydrogenated models.The electronic densities of states (DOS) for the pristine and hydrogenated a-P models (Figure 3b) were computed using hybrid DFT.All three pristine models show mid-gap defect states between the valence-band maximum and the conduction-band minimum.The band-resolved charge densities for the highest occupied bands (Figure 3c) show that these mid-gap states are mostly caused by non-bonded electrons around the twofold-connected atoms and by the defect states at the fourfold-connected atoms in pristine a-P models.Upon hydrogenation, all defect states disappeared, increasing the band gaps in all three hydrogenated models compared to their pristine counterparts.The 3H model has a slightly larger predicted gap (2.06 eV) than 1H (2.01 eV) and 2H (1.86 eV), suggesting that structural ordering might play a role in opening up the electronic band gap in a-P.The computed band gaps for our a-P models agree well with various, previously reported experimental results for amorphous red phosphorus (1.42 to 2.07 eV). [38]panding upon prior studies of crystalline allotropes, [5b, 39] we also calculated the electronic DOS, at the hybrid-DFT level, for all crystalline phases discussed in this work (Figure S3)allowing for side-by-side comparison with our a-P models.A wide range of band gaps was found for the various phosphorus modifications: white P has the largest gap (3.80 eV) among the crystalline allotropes, consistent with its molecular nature, whereas only a small gap exists in black P (0.23 eV), in agreement with experimental data (≈ 0.3 eV). [15,40] y contrast, violet and fibrous phosphorus have moderate gap sizes (2.26 and 2.30 eV, respectively), close to that of 3H, which shows a relatively more ordered network than the other two hydrogenated models.
The band gap sizes of the nanorod structures are in the same range as for other forms of P containing cage-like motifs (viz.2.01-2.26eV, Figure S3).
The Crystal Orbital Hamilton Population (COHP) technique allows one to quantify the bond strength and understand the bonding features, based on orbital interactions from the hybrid DFT calculations.We calculated -COHP curves for all P-P bonds in the crystalline and amorphous modifications mentioned above; negative COHP values stand for stabilisation (bonding).The integrated -COHP (referred to as -ICOHP) of a given bond, up to the Fermi level, EF, has been used to quantify the strength of chemical bonds in similar systems. [41]Figure 4a shows the relative bond strength, normalised to the shortest bond in crystalline (violet) phosphorus, for all crystalline modifications and up to 2.7 Å.An empirical measure, introduced by Pauling in the 1940s, can be used to describe this relation between bond length and strength: [42] 1 −  =  log 10  in which 1 is the single P-P bond length based on tabulated covalent radii (2.22 Å), [42] and  is the length of different P-P bonds in the system. is the normalised -ICOHP value, and  is a fitting coefficient.During fitting, data for white phosphorus were not included, since its bonding nature, involving P4 molecules, clearly deviates from that of the other crystalline modifications.Bonds longer than 2.7 Å were also ignored, as Pauling's formula is only expected to describe well the region of relatively strong covalency.
We found that the fitted Pauling relation for black, violet, and fibrous phosphorus almost overlaps with the one fitted for P nanorods, suggesting similar bonding in all these crystalline modifications.The same analysis was performed for all P-P bonds in the three pristine a-P models (Figure 4b).The -ICOHP data for a-P scatter in a wide range due to more complex bonding environments; the bonds in the single P4 unit in model 1 are located away from the others, in the same region as found for crystalline white phosphorus.Fitting Pauling's relation based on the rest of the bonds suggests slightly different bonding in a-P than in the crystalline phases: the bond strength diminishes more quickly with bond length in a-P, perhaps indicating slightly lower chemical stability of the extended network.
Energy-resolved -COHP plots provide more detailed "fingerprints" of bonding, as shown in Figure 4c.Despite stable bonds being formed in white phosphorus (within the P4 tetrahedra), evident from stabilising areas below EF (-COHP > 0), the overall bond strength as measured by -ICOHP is lower than in the other forms.By contrast, a small antibonding interaction below EF is observed for black phosphorus (-COHP < 0).Such interactions were reported for P-P bonds in structurally related compounds, viz.(Li-intercalated) phosphorene [43] and BaP4Te2, [44] and in various chalcogenides, such as the iso-valent-electronic GeTe [45] and related Ge-Sb-Te alloys. [46]Hence, antibonding areas below EF do not necessarily suggest poor stability, especially as there are no antibonding interactions directly at EF, and the -ICOHP value is large.
Violet and fibrous phosphorus show only bonding areas below EF, and so does n1 (with P8 building units similar to those in the violet and fibrous forms); more complex cage motifs (e.g., P10 and P12) might result in marginal antibonding regions at the valence-band edge in n2 and n3.Similar interactions are seen for a-P: there are small antibonding regions not only for the defect states in the gap, but also below the top of valence band.The former is largely attributed to mid-gap defect states; the latter may be due to complex cage motifs (such as in n2 and n3) that are less "ideal" than the P8 / P9 building units in the violet and fibrous forms.

Conclusions
We have created structural models of amorphous red phosphorus by combining machine-learning-driven molecular dynamics and first-principles computations.Our a-P models are energetically intermediate between white phosphorus and the other crystalline modifications.We find that the details of their energetic stability depend on the degree of structural ordering: the less stable amorphous model 1 is rather disordered and resembles a random network, whereas the more stable model 3 is structurally similar to crystalline tubular phases (e.g., violet and fibrous phosphorus).Our work completes the first-principles investigation of the stability range of phosphorus modifications, [21] having added the challenging case of a-P to the picture.We also created models of hydrogenated models a-P, thereby revealing the impact of defect states on the electronic properties of pristine a-P.Our chemical-bonding analyses quantified the relation between bond length and strength, indicating slightly weaker covalent bonding in a-P than in its crystalline counterparts.Whilst the present work has focused on the fundamental structural chemistry of a-P, the structural models provided herein are amenable to further DFT investigations: as one example, we envision first-principles computational studies of Na-ion insertion in chemically more complex a-P-based battery materials.
placed perpendicularly to the two chemical bonds of that P atom, on the side with more open space (i.e. with the lowest sum of distances to all neighbouring P atoms in the local environment), forming a short P-H bond.The initial bond length was set to 1.4 Å, a typical value found in phosphorus hydrides. [54]The resultant structures were then computationally optimised using PBE+D3.
DFT computations.Structural relaxations and single-point calculations were performed using the Vienna Ab initio Simulation Package (VASP) [55] with projector augmented-wave (PAW) [56] pseudopotentials.The PBE+D3 method [24,57] was used for structural relaxations, whereas the HSE06 hybrid functional [58] with many-body dispersion (MBD) corrections [25] was used in the subsequent static computations, based on pre-converged PBE wave functions.The plane-wave energy cut-off was 500 eV, the energy tolerance for SCF convergence was 10 -7 eV per cell, and the force tolerance for structural relaxation was 0.01 eV Å -1 .Gaussian smearing with a width of 0.05 eV was used to determine partial occupancies during SCF cycles.In structural relaxation, all lattice parameters and atomic coordinates were optimised.A k-point grid with a largest allowed spacing of 0.1 Å -1 along each reciprocal lattice vector was used for structural relaxation for crystalline modifications, and a spacing of 0.15 Å -1 (i.e., a lower k-point density) was used for the following HSE06 computations of the electronic structure and bonding, to reduce computational cost.We verified that changing between both settings did not change the predicted band gap of black phosphorus (0.23 eV) within the quoted accuracy.For the a-P and a-P:H models, all data reported are from Γ-point-only computations.
Chemical-bonding.The Local-Orbital Basis Suite Towards Electronic-Structure Reconstruction (LOBSTER) code [59] was used to project the self-consistent wave functions onto an auxiliary basis of localised, atom-centered orbitals (3s and 3p on each P atom), enabling crystal orbital Hamilton population (COHP) analysis. [60]sualisation.The structural models in Figure 1a and 3a, as well as the fragments shown in Figure 2a, were visualised using OVITO; [61] the visualisations of local motifs and band-resolved charge density in Figure 3c were created using VESTA. [62]

Figure 1 .
Figure 1.Structural models of amorphous phosphorus.(a) Structures of a-P as obtained from ML-driven simulations and subsequent DFT structural optimisation.We generated three independent structural models containing 248 atoms each, using the same simulation protocol in each case.The structures are labelled 1 to 3, with an increasing degree of structural order.Colour-coding indicates the coordination numbers, viz.N = 2 (blue), N = 3 (pale blue), and N = 4 (pink), determined using a 2.4 Å cut-off.(b) Relative abundance (per-atom count) of cluster fragments in the a-P models (circles), as compared to the result of a similar but larger-scale 1,984-atom simulation in Ref. [18] (squares).Solid and dashed lines are guides to the eye.Sketches of selected fragments are shown, following Ref.[20].

Figure 2 .
Figure 2. The energetic stability range of phosphorus modifications, including the well-known white, black, violet, and fibrous forms, as well as structural models of Pfitzner's nanorods (n1 to n3) based on Refs.[7] and [21], and the new a-P structural models 1 to 3 generated in the present work.(a) Selected, local structural fragments in our a-P models compared to the building units of nanorods.(b) Computed energies, given relative to the standard state, viz.white phosphorus, and based on DFT computations (HSE06+MBD).

Figure 3 .
Figure 3. Hydrogen in a-P and its effect on the electronic structure.(a) Structural models for a-P:H, 1H to 3H.To generate these structural models, four-fold-connected atoms were removed and resulting two-fold-connected atoms saturated with hydrogen atoms, followed by relaxation at the PBE+D3 level.Consequently, these structures are fully connected according to standard valence rules.(b) Computed electronic density of states (DOS) plots for pristine and hydrogenated a-P, using HSE06.Dashed lines indicate the DFT-computed Fermi level, EF, here taken to be located at the top of the valence band in the hydrogenated models.The energy scales (y-axis) for the pristine models were each shifted such that the zero energy represents EF of the corresponding hydrogenated model.(c) Visualisation of the band-resolved charge density for the highest occupied band (grey isosurfaces) of a defect state in model 1 (upper panel) and its corresponding state (lower panel) in model 1H after hydrogenation.The coordination defect (N = 2) in 1 and the same atom (now fully connected; N = 3) in 1H, together with the added hydrogen atom, are marked by arrows.Except for the marked atoms, all other atoms are fully connected (N = 3) in the structural fragments shown.Different isovalues are used (viz.0.006 in 1; and 0.001 in 1H), to only highlight the charge distribution around the marked atoms.Despite the lower isovalue used for 1H, no localised electrons were found on visual inspection, corroborating the removal of defect states upon hydrogenation.

Figure 4 .
Figure 4.Chemical bonding in crystalline and amorphous phosphorus.Bond-length-bondstrength correlations in the different modifications are obtained from integrated crystal orbital Hamilton population (ICOHP) analysis of: (a) crystalline modifications, including the nanorod models n1-n3; and (b) a-P models.In these plots, ICOHP values are all normalised using a factor of 6.07 eV, such that the shortest bond (2.18 Å) in crystalline violet phosphorus has a value of 1.0.White phosphorus stands apart from the other crystalline models, whereas the rest of the crystalline structures and nanorods all follow a similar trend.Curves representing Pauling's relation are shown by dashed lines (see text).We note that the bonds in the single P4 unit in the a-P model 1 are located in the same region as those for crystalline white P in panel (a); in both cases, these data points have not been included in the fits.(c-e) Energy-resolved -COHP plots, emphasising the close similarity between a-P, nanorods, and violet and fibrous phosphorus.The positive and negative signs on the abscissae indicate the bonding and antibonding regions, respectively.The dashed horizontal line indicates the Fermi level, EF.

Figure S1 .
Figure S1.Generating a-P structural models with machine-learning-driven simulations.(a) The temperature profile of the melt-quench process used to generate the amorphous models in this work (as described in the Methods section): a simulated disordered liquid phase is cooled over about 10 ns, or 10 million simulation steps.(b) Counts of coordination numbers during the melt-quench process.In the liquid phase, most atoms are two-(≈ 49 ± 4%) and three-fold (≈ 29 ± 6%) connected.During quenching, a network of mostly three-fold-connected atoms emerges, and in the final structures (indicated by markers), almost all of the atoms have N = 3 neighbours.(c-d) Calculated radial distribution function (RDF, panel c) and angular distribution function (ADF, panel d) for the three pristine models of a-P (1 to 3) generated in the present work, compared to those calculated from the larger-scale model in previous work [Y.Zhou, W. Kirkpatrick, V. L. Deringer, Adv.Mater.2022, 34, 2107515].