Open Access Article
Laura Bonometti
a,
Giuseppe Sansonea,
Marcos Rivera-Almazo
a,
Denis Usvyatb,
Antti J. Karttunen
*c and
Lorenzo Maschio
*d
aDipartimento di Chimica, Università di Torino, Via P. Giuria 5, 10125 Torino, Italy
bInstitut für Chemie, Humboldt-Universität zu Berlin, Brook-Taylor-Str. 2, D-12489 Berlin, Germany
cDepartment of Chemistry and Materials Science, Aalto University, Kemistintie 1, 02150 Espoo, Finland. E-mail: antti.karttunen@aalto.fi
dDipartimento di Chimica, NIS Centre, Università di Torino, Via P. Giuria 5, 10125 Torino, Italy. E-mail: lorenzo.maschio@unito.it
First published on 6th November 2025
In this work, we investigate through quantum–mechanical calculations the relative stability of white-γ, white-β, fibrous red, violet, and orthorhombic black phosphorus allotropes, a longstanding yet challenging problem. We performed DFT-D3 calculations with the CRYSTAL code as well as periodic local second-order Møller–Plesset perturbation theory (p-LMP2) calculations. DFT and Spin-Component Scaled p-LMP2 place violet phosphorus as the thermodynamically most stable allotrope both at 0 K and at 298 K, yet within a tiny margin from the black phosphorus. Pure p-LMP2 suggests that black phosphorus is the most stable, although its accuracy may be affected by the narrow band gap in this material.
The white phosphorus modifications are based on P4 tetrahedra, and there are two different crystalline forms, namely β and γ.10–15
The red phosphorus class exhibits the greatest structural diversity2 and only two crystalline allotropes are known: fibrous red (red-IV) and violet (Hittorf's, red-V) phosphorus.16–23 They are made of tubular P units that form complex double layers by connecting perpendicularly in violet phosphorus and in parallel in fibrous red phosphorus. The layers are then stacked in the three-dimensional crystal structure, stabilized by weak van der Waals forces.
Several modifications also exist in the black phosphorus class. At ambient temperature and high pressures (20 °C; 12 kbar), crystalline orthorhombic black phosphorus can be obtained, which shows a layered structure composed of puckered six-membered P rings. The layers are kept together also by van der Waals interactions.3,4,11,24
Bridgman's discovery of the orthorhombic black phosphorus with layered structure encouraged further studies to determine the structural and thermodynamical relations of the phosphorus allotropes.11 Orthorhombic black phosphorus25,26 is typically described as the thermodynamically most stable form of the element at NTP conditions.3,4 However, re-tracing the main historic studies on the topic, Aykol et al.27 have suggested that there is little direct evidence to support this claim.11,19,24,28–35 Nilges et al.1 have discussed the difficulties in obtaining hard thermodynamic evidence on the stability order of bulk phosphorus allotropes from temperature-dependent vapor-pressure measurements.
The structural variety of phosphorus allotropes and the experimental difficulties in pinpointing the thermodynamic stability order have led in several systematic computational investigations within this field with density functional theory (DFT).27,36–39 New phosphorus modifications are also being proposed and predicted theoretically, confirming the lively interest in this element.40 Furthermore, machine-learning interatomic potentials based on DFT plus many-body dispersion approach have enabled the investigation of both crystalline and amorphous phosphorus modifications in size scale that was inaccessible few years ago.41–44
From a computational point of view, one of the most challenging aspects of the thermodynamical comparison of various phosphorus allotropes is the coexistence of covalent bonding and weak dispersive interactions. Standard DFT methods do not describe the weak van der Waals interactions in phosphorus allotropes properly and dispersion-corrected DFT methods are needed.45–50 While several dispersion-corrected DFT approaches are available, they need to be benchmarked carefully to ensure that all classes of phosphorus allotropes are described with equal accuracy.
In this work (Fig. 1), we investigate the thermodynamics of bulk crystalline phosphorus allotropes that exist at ambient pressure, deliberately excluding the high-pressure phases, such as those of black phosphorus,36,51 which lie beyond the scope of this study. The studied phosphorus allotropes are listed in Table 1 and their crystal structures are illustrated in Fig. 2. Dispersion-corrected DFT methods are used to fully optimize the structures and determine the energetic and thermodynamic stability of the phosphorus allotropes, together with other thermodynamic quantities such as the entropy and heat capacity. For an ab initio description of dispersion we also utilize periodic local second-order Møller–Plesset perturbation theory (p-LMP2).
![]() | ||
| Fig. 2 Structures of the phosphorus allotropes listed in Table 1. | ||
We believe that a coherent dataset of structures, energies and thermodynamic functions provides a robust benchmarking reference, which complements previous theoretical efforts to determine the most stable phosphorus allotrope.
For sampling the reciprocal space, Monkhorst–Pack type k-meshes were used: 4 × 4 × 6 for mS8, 6 × 2 × 2 for aP24, 3 × 3 × 5 for aP42, 4 × 4 × 1 for mP84 and 6 × 6 × 6 for oS8.59 Both nuclear coordinates and lattice parameters were relaxed for the geometry optimization of all studied structures, conserving the space group symmetry. Default DFT integration grids and optimization criteria of CRYSTAL23 were applied. The five thresholds Ti, which control the truncation criteria of the Coulomb and exchange infinite lattice series, have been set to 8 (T1–T4) and 16 (T5). To properly describe the thermodynamic properties of the allotropes, we sampled phonon modes outside the Γ-point by using phonon supercells as implemented in the CRYSTAL code.54,60 The applied phonon supercells and the k-meshes used in the phonon calculations are reported in SI.
Then, we considered the optimized geometries at PBE0-D3(ZD)/TZVPP level of all the studied systems for single-point calculations with the periodic local second-order Møller–Plesset perturbation theory (p-LMP2) and its spin-component-scaled variant61 (p-SCS-LMP2).62–65 As the HF reference for p-LMP2, we used the AO-based periodic HF of the CRYSTAL code.54
In the p-LMP2 method the occupied manifold is represented by local Wannier functions (WFs),66,67 allowing one to exploit the short-range nature of electron correlation and to achieve a nearly linear scaling with the number of atoms per cell. As virtual functions p-LMP2 employs the so-called orbital specific virtuals (OSVs),68,69 which provide an efficient and compact representation of the pair-specific virtual space. The two-electron integrals in p-LMP2 are approximated via local density fitting.70
| Param. | mS8 | aP24 | aP42 | mP84 | oS8 |
|---|---|---|---|---|---|
| a | 9.25 (+0.9%) | 5.48 (0.0%) | 12.27 (+0.6%) | 9.22 (+0.1%) | 3.30 (−0.1%) |
| b | 8.25 (−1.1%) | 10.86 (+0.6%) | 13.00 (+0.1%) | 9.13 (−0.2%) | 10.61 (+1.3%) |
| c | 5.45 (+0.4%) | 11.08 (+1.1%) | 7.11 (+0.5%) | 22.65 (+0.2%) | 4.42 (+0.9%) |
| α | — | 93.9 (−0.4%) | 116.9 (−0.1%) | — | — |
| β | 90.5 (+0.2%) | 99.7 (0.0%) | 106.4 (+0.1%) | 106.0(−0.1%) | — |
| γ | — | 101.2 (+0.5%) | 97.9 (0.0%) | — | — |
| V | 415.8(+0.1%) | 634.4 (+1.7%) | 922.8 (+1.3%) | 1832.9 (+0.2%) | 154.8 (+1.9%) |
The temperature-dependent Gibbs free energy curves at ambient pressure for the solid-state phosphorus allotropes and for the isolated (gas-phase) P4 molecule are shown in Fig. 3, where also experimental phase transitions are reported to be compared with our results. Some of the experimental results could not be simulated as they involve amorphous phosphorus allotropes such as as amorphous red phosphorus (red-I) or α-white modification.
![]() | ||
| Fig. 3 Gibbs free energy curves of the studied phosphorus allotropes and the gas-phase P4 molecule. The zero level is the Gibbs free energy of P4 at 800 K. The symbols highlight the calculated crossovers between the Gibbs free energies of the allotropes. Some experimental phase transitions are reported as vertical lines, according to the ref. 1 for mP84 and oS8 allotropes transitions to P4(g); to ref. 12 for experimental transition from γ to β white modification; to ref. 34 for the transition from aP42 to P4(g); to ref. 3 for transition from red-I to P4(l) and to ref. 4 for all the remaining transitions. The inset highlighted by the green rectangle is shown in high detail in Fig. 4. Red-I refers to the amorphous red phosphorus. | ||
First we note that, as concerns the white modifications mS8 and aP24, the calculated results agree reasonably well with the experiment. The computed values indicate that these allotropes are substantially less thermodynamically stable than the black or red/violet phosphorus, which is not surprising as white phosphorus modifications are weakly bound molecular crystals. Experimentally this manifests itself in relatively low temperatures of the phase transition from the aP24 to the amorphous red-I phase and later from red-I to the liquid and gas phases. Since the amorphous phases are not directly accessible computationally, the calculation considers a direct transition from aP24 to the gas phase. The predicted temperature of such a transition is slightly underestimated compared to the experimental temperatures, which may suggest that the PBE0-D3 underestimates the dispersive interactions in white phosphorus (see also discussion below on the p-LMP2 results) or that the harmonic approximations fails for the soft intermolecular modes (or both).
According to our calculations the two allotropes of white phosphorus are energetically very close to each other at all temperatures, within sub-kJ mol−1 per atom, see Table 3. At high temperatures the aP24 modification is slightly more stable than mS8, at low temperatures the stability order reverts. This trend agrees with the experimental phase transition from mS8 to aP24 by heating. Notably, the calculated phase-transition temperatures align fairly well with experiment (see Fig. 3), which is rather remarkable given the underlying harmonic approximation and the relative crudeness of the D3 dispersion model.
| Allotrope | ΔE | ΔG0 | ΔG298 | ΔG500 |
|---|---|---|---|---|
| mS8 | 20.0 | 19.2 | 16.4 | 14.0 |
| aP24 | 20.3 | 19.4 | 16.3 | 13.6 |
| aP42 | 0.4 | 0.4 | 0.1 | −0.1 |
| mP84 | −0.1 | −0.1 | −0.3 | −0.5 |
| oS8 | 0.0 | 0.0 | 0.0 | 0.0 |
Next we focus on the other three allotropes and the temperature region of the phase transitions to the gas phase, which is shown in detail see Fig. 4. Experimentally the temperatures of such phase transitions lie in a very small window, with the temperature for the black phosphorus' transition being the highest. The calculations do not exactly reproduce this order, most likely due to the harmonic approximation, possibly in combination with the approximations within DFT and the dispersion model employed. Indeed, given the presence of soft and likely anharmonic inter-layer or inter-fiber modes, it is not surprising that the harmonic approximation does not correctly resolve the tiny differences in the Gibbs free energies at high temperatures. Nevertheless, and more importantly, our calculations reproduce the proximity of the phase transition temperatures and furthermore predict the temperature range for these transitions very accurately. In Fig. 3 and Table 4 all relevant experimental and calculated phase transitions temperatures are illustrated and listed.
![]() | ||
| Fig. 4 Zoomed-in portion of Gibbs free energy curves in between 650 and 800 K of os8, aP42 and mP84 allotropes and the gas-phase P4 molecule. The zero level is the Gibbs free energy of P4 at 800 K. The symbols highlight the calculated crossovers between the Gibbs free energies of the allotropes. Experimental phase transitions are reported as vertical lines, according to the ref. 1 for mP84 and oS8 allotropes transitions to P4(g) and to ref. 34 for transition from from aP42 to P4(g). | ||
For a further insight in the accuracy of our calculated Gibbs free energies, we compared our results with the experimental heat capacity and entropy data. Stephenson et al.19 reported the constant pressure heat capacities of phosphorus allotropes from about 15 up to 300 K, together with entropies at low temperatures and at 298 K. Other data are provided by ref. 80 and 81. Paukov et al.82 published data for the oS8 allotrope, but according to Stephenson et al. their data differ only about 1%.19 Comparisons of the entropies and heat capacities are reported in Tables 5 and 6 respectively. The entropy data show a good overall agreement between the calculations and the experiments. The calculated trends of the heat capacities also agree well with the experimental ones, see Fig. 5. In particular, the low-temperature behavior of the mS8 and aP24 allotropes is well reproduced. Above 158 K, mS8 converts to aP24, which in turn transforms to amorphous white phosphorus at about 196 K.12 Right in this region, from about 150 K, the slopes of the experimental curves change in a way that indicate the upcoming phase transition.
![]() | ||
| Fig. 5 Left: heat capacities at constant volume predicted in this work. Right: experimental heat capacities at constant pressure.19 Three different forms of white phosphorus, white-β, white-γ, and amorphous white are present in the temperature range, but they are not distinguished in the experimental data. | ||
Finally, we look at the low-temperature regime for the black and red/violet phosphorus allotropes. As the temperatures drops, the oS8 is predicted to become more stable than aP42. However it still remains slightly above mP84 at any temperature (see Fig. 3 and Table 3), with a minute energy difference of −0.1 kJ mol−1 per atom in favor of the latter at T = 0. This result is in line with the predictions of other theoretical investigations. Generally the sign of the relative stability between oS8 and mP84 may vary depending on the functional employed.27,38 However calculations using high-level dispersion models, such as the random phase approximation (RPA)27 or HSE06+MBD (many-body-dispersion correction)44 both predict mP84 to be more energetically stable, albeit by a tiny amount.
| Allotrope | ΔEp-LMP2 | ΔEp-SCS-LMP2 |
|---|---|---|
| mS8 | 12.74 | 9.78 |
| aP24 | 13.81 | 10.38 |
| aP42 | 3.61 | 1.99 |
| mP84 | 1.49 | −0.19 |
| oS8 | 0.00 | 0.00 |
Interestingly, p-LMP2 predicts the black phosphorus to be the most stable allotrope. However the p-LMP2 result should be interpreted with caution as MP2 is known to overestimate dispersion in polarizable systems in general. Furthermore, MP2 is not able to capture Coulomb screening, which effectively weakens the long-range part of the dispersive interactions in bulk materials. Such screening effects are strong in narrow gap semiconductors, which is an additional source for overestimation of the interaction energy by MP2. Black phosphorus is a particular example, where pure p-LMP2 is known to significantly overestimate the exfoliation energy.49 We note that finite-temperature effects on electronic band gaps arise from electron–phonon coupling and can modulate correlation energies. Explicit inclusion of electron–phonon coupling for bulk black phosphorus is, however, beyond the computational scope of the present study. Nonetheless, experimental data indicate a temperature dependence of the band gap of about
meV T−1 in between 100 and 300 K83 and such a small absolute variation is unlikely to qualitatively affect our conclusions regarding MP2 overbinding. The violet phosphorus still being a narrow gap semiconductor, has a wider gap than a black phosphorus (e.g. the PBE0 gaps in our calculations: 0.90 eV for oS8 and 2.77 eV for mP84). Therefore p-LMP2 is expected to be less inaccurate for the violet phosphorus than for the black one, resulting in an artificial relative overstabilization of the latter.
The SCS-MP2 technique was introduced84 to repair the MP2 overbinding in highly polarizable systems. Therefore, p-SCS-LMP2 results for our systems are expected to be more balanced. And indeed the relative stability with p-SCS-LMP2 is again reversed, making the violet-phosphorus the most stable allotrope (yet again by a just −0.2 kJ mol−1 per atom). That is particularly interesting, as p-SCS-MP2 does not capture the Coulomb screening effects either50 and is expected to still be somewhat biased towards black phosphorus, although not as much as p-LMP2. And yet, it already favors the violet phosphorus, in line with the other high-level methodologies like RPA27 or HSE06 + MBD44 or HSE06 − MBD,44 which do include the screening effects.
As concerns other studied allotropes, the p-SCS-LMP2 results qualitatively follow the PBE0-D3 order of relatives stabilities. There is however a quantitative disagreement for the stability of the white phosphorus modifications: p-SCS-LMP2 predicts the white phosphorus to be energetically much closer to the black phosphorus than PBE0-D3. This result may be another indication that PBE0-D3 noticeably underestimates binding in white phosphorus, as suggested by the discrepancy in the predicted temperature of the phase transition to the gas phase (see above). Other theoretical estimates of the energy difference between the white and black phosphorus are about 15 kJ mol−1 per atom,1,27,38 which is somewhat in between the p-SCS-LMP2 and PBE0-D3 predictions.
DFT and p-SCS-LMP2 at 0 K predict that the violet allotrope is energetically more stable than black phosphorus, but only minutely: by just 0.1–0.2 kJ mol−1 per atom. The p-LMP2 prediction reverts this order, but most likely due to the gross overestimation of dispersion in black phosphorus due to its very narrow gap.
Even if normalizing the energy difference between violet and black phosphorus by a characteristic unit like P4 or a unit cell of black phosphorus (eight atoms), it still remains clearly below the error margin of the methods employed here or in the other studies. Therefore, with confidence one can only state that these two allotropes are nearly isoenergetic. However, given the agreement between very different and relatively high-level theoretical approaches, one can cautiously speculate that it is the violet phosphorus being the most stable allotrope at low temperatures. For the other studied modifications the predicted order – fibrous red → γ-white → β-white—agrees with experimental data and other calculations.
At higher temperatures, PBE0-D3 with harmonic approximation reasonably accurately captures the temperature of the phase transition between the γ-white and β-white phosphorus, but somewhat underestimates their hypothetical sublimation temperatures. This suggests that PBE0-D3 may noticeably underestimate the formation energy of the molecular crystals of white phosphorus. For the black and red/violet phosphorus the predictions of their sublimation temperatures, which are very close to each other, are quite accurate. However, the order of the temperatures is not reproduced probably due to the deficiencies of the harmonic approximation.
| This journal is © The Royal Society of Chemistry 2025 |