How to evaluate aromaticity under pressure? Benzene as a benchmark system

Abstract

Aromaticity is a fundamental concept in chemistry, central to understanding chemical stability, reactivity, and electronic structure. While traditionally examined under ambient conditions, recent studies have begun to probe its behavior in unconventional regimes, reflecting the diversification of reaction environments. Here, we investigate aromaticity under hydrostatic pressure, extending the range of external perturbations considered so far in aromaticity research. Using state-of-the-art quantum chemical methodologies designed to simulate pressure at the single-molecule level, we systematically analyze the evolution of aromaticity in benzene, the prototypical aromatic compound. Structural and electronic indices consistently indicate a modest loss of aromaticity under pressure, while magnetic descriptors point towards an enhanced aromaticity under pressure. This underscores the enduring need for a multidimensional framework consisting of complementary descriptors. During the analysis, particular emphasis is placed on the performance, adaptability, and limitations of widely used aromaticity indices under compression, leading to proposed refinements and guidelines for their reliable application. The results are interpreted within different theoretical perspectives on aromaticity, providing new insights into the robustness of aromaticity descriptors under extreme conditions and outlining avenues for future research.

---

Introduction

Aromaticity is a cornerstone concept in chemistry, underpinning our understanding of stability, reactivity, and electronic structure across organic, inorganic, and biological systems.^1–3 Despite ongoing debates regarding its definition,^4,5 it remains a powerful framework for rationalizing chemical behavior and guiding the design and discovery of novel molecules and materials.⁶ Its enduring relevance is reflected in the growing bibliometric impact,¹ with over a thousand publications on aromaticity in 2024 alone.

Aromaticity, however, is not a direct quantum observable, nor are most of its associated electronic structure consequences. Similar to other fundamental concepts such as chemical bonding or oxidation state, it must be defined by convention,⁷ and numerous methods have been developed to quantify aromaticity based on energetic,⁸ geometric,⁹ electronic¹⁰ and magnetic criteria.¹¹ Each descriptor provides only an indirect assessment of a particular facet of aromaticity, and none is universal or without limitations.¹² Accordingly, it is advisable to employ a set of complementary descriptors grounded in distinct criteria.^3,13 When multiple indices are considered, conclusions are generally more reliable, particularly when all criteria converge on the same trend.¹⁴ Nonetheless, the multidimensional character of aromaticity should not be invoked as a justification for accepting any index as valid, regardless of the results obtained.^15,16

Since its inception, benzene has served as the paradigm of aromaticity. Its distinctive features, remarkable thermodynamic stability, cyclic electronic delocalization, bond-length equalization, characteristic magnetic responses and unusual reactivity, naturally led to their association with the concept of aromaticity.¹⁷ For benzene, all descriptors consistently assign a high degree of aromatic character, and most indices, though not all, identify it as the most aromatic species. Yet discrepancies emerge when comparing other systems, as correlations among indices depend strongly on both the compound set and the aromaticity indices considered.^3,18–20 Benzene therefore provides a natural benchmark, but a deeper understanding of aromaticity requires testing the descriptors under external perturbations beyond ambient conditions.^21,22

Recent studies have shown that electric fields can modulate the aromatic character of molecules, even suggesting tunability via the field intensity.^23,24 Additionally, magnetic fields have been shown to be able to deform the benzene molecule, although significant effects are only seen at very high field strengths.^25,26 Building on this perspective, we extend the investigation to another fundamental perturbation: pressure. In this work, we systematically explore how the aromaticity of archetypal benzene evolves under compression, employing a diverse set of descriptors grounded in distinct criteria. To the best of our knowledge, very few studies have examined how aromaticity varies under compression, and those studies primarily focus on molecular crystals.^27–29 Yet, pressure often challenges our understanding of fundamental chemical concepts,^30–32 and in this context, we address a central question: how can aromaticity be reliably evaluated under pressure?

To address this question, we employ state-of-the-art computational approaches to simulate pressure at the single-molecule level.³³ While high-pressure experiments necessarily involve molecules embedded within a chemical environment, modeling the response of isolated molecules allows intra- and intermolecular effects to be disentangled, thereby offering fundamental insights into molecular behavior under compression.

In this study, we employ two complementary methods to investigate the pressure-induced properties of benzene. First, we apply the eXtended Hydrostatic Compression Force Field (X-HCFF) approach, which imposes compressive forces perpendicular to the discretized van der Waals surface of a molecule by manipulating the nuclear gradient during geometry optimization.^34,35 This procedure provides direct access to optimized molecular geometries at a user-defined pressure, enabling efficient simulation of pressure-induced structural changes.

Second, we employ the Gaussians On Surface Tesserae Simulate HYdrostatic Pressure (GOSTSHYP) approach.^36,37 In addition to shrinking the nuclear framework, GOSTSHYP directly compresses the electron density through a smooth field of Gaussian potentials distributed over a suitably selected molecular surface.³⁸ Although more demanding both theoretically and computationally than X-HCFF, GOSTSHYP yields more physically accurate electronic structures at elevated pressures.³⁶ Applying both X-HCFF and GOSTSHYP therefore enables us to disentangle nuclear and electronic contributions to pressure-induced changes in aromaticity.

To complement these single-molecule simulations, we also perform high pressure simulations on a benzene crystal.^39,40 In this extended model, a benzene molecule is embedded within the anisotropic environment created by neighboring benzene molecules, offering a more realistic representation of experimental conditions. This approach contrasts with the isotropic pressure applied in X-HCFF and GOSTSHYP, providing complementary insights into how pressures in the GPa regime influence the aromaticity of benzene.

The paper is structured as follows. First, the computational methodology is briefly described. Next, the behavior of selected aromaticity descriptors under pressure is examined, with the indices grouped by category. Each group is introduced with a concise description of the relevant indicators, followed by an analysis of the aromaticity trends. Where applicable, adaptations to the indices and the implications of pressure are discussed. Finally, the extent to which the different criteria converge in describing pressure-induced changes in aromaticity is evaluated.

Computational details

To evaluate the aromaticity of benzene under pressure, geometry optimizations were first performed using the X-HCFF and GOSTSHYP methods implemented in the Q-Chem 6.0 software.⁴¹ All calculations are based on Kohn–Sham Density Functional Theory (KS-DFT), employing the B3LYP functional⁴² with D3 dispersion corrections⁴³ and the cc-pVTZ⁴⁴ basis set. Diffuse functions were deliberately omitted to ensure that, during compression with GOSTSHYP, no electron density extends beyond the defined cavity.

In GOSTSHYP, the molecular surface plays a critical role in compressing the electron density, directly influencing the aromaticity descriptors.⁴⁵ To approximate the molecular shape under ambient conditions, the 0.001 e bohr⁻³ isodensity surface was used. Based on the van der Waals radii of carbon and hydrogen, the cavity was constructed as a point set along an angle-dependent grid centered around each atom. The quality of the resulting molecular surface can then be quantified as the root-mean-square error RMSE(C) of the minimum distance from a point to the isodensity surface.

here, C denotes the cavity, and i runs over the discretization of the isodensity surface. and correspond to the point positions along an angular grid of the isodensity surface and the cavity, respectively. In the next step, the deviation from the desired shape, as defined by the isodensity surface, was minimized. This was achieved through an iterative minimization of RMSE(C) using the steepest descent method.⁴⁵ Both the calculation of the RMSE relative to the isodensity surface and the associated steepest descent algorithm were implemented in a Python 3.11 script.⁴⁶ The evolution of the RMSE(C) and the van der Waals radii for carbon and hydrogen during the minimization are shown in Fig. S1 along with optimized radii values. In Fig. S2, the 0.001 e bohr⁻³ isodensity surface and the optimized cavity are visualized.

To exclude the possibility that long-range interactions influence the aromaticity trends under pressure,^47,48 additional pressure calculations were performed on the gas-phase molecule with GOSTSHYP using the range-separated CAM-B3LYP functional⁴⁹ and the optimized van der Waals radii, while keeping other settings unchanged. On this basis, structural indices, ring critical point (RCP) properties and electronic descriptors were evaluated, all of which reproduced identical pressure-dependent aromaticity trends as obtained with B3LYP (Table S6 and Fig. S37). Furthermore, a technical limitation of the Q-Chem code prevents the calculation of NICS indices and magnetic descriptors using the CAM-B3LYP functional. Consequently, B3LYP + D3 was employed to ensure internal consistency across all calculated indices.

To compare the aromaticity of a benzene molecule in the gas phase with that in a crystalline environment, plane-wave-based DFT calculations under pressure were performed on the crystal structure using the Vienna Ab initio Simulation Package (VASP) 3.6.0.^50,51 To describe the electron interactions, the projector augmented-wave (PAW) method^52,53 was employed together with the GGA-PBE exchange-correlation functional,⁵⁴ using appropriate PAW pseudo potentials for carbon and hydrogen. The wavefunctions were represented using a plane-wave energy cutoff of 400 eV, while the Brillouin zone was sampled using a Gamma-centered 3 × 3 × 2 k-point mesh. After crystal optimization, the geometry of a single benzene molecule was extracted at each pressure for further analysis.

Using the above pressure models, five main approaches were employed to perform pressure calculations prior to further analysis. The “Gas – X-HCFF” and “Gas – GOSTSHYP” methods apply the standard X-HCFF and GOSTSHYP procedures, respectively, in both the geometry optimization and electronic structure calculations. The “Crystal Geo” method extracts the geometry from the crystal structure and computes the corresponding gas-phase electronic structure. To account for the electron cloud compression, the GOSTSHYP method can be used in the final single-point calculation on this geometry, referred to as “Crystal Geo – GOSTSHYP”. Finally, the “Gas – GOSTSHYP Geo” method was introduced to disentangle the effect of geometry distortion and electron-cloud compression in GOSTSHYP by performing a gas-phase single-point calculation on the optimized geometry obtained from the GOSTSHYP simulation.

To appropriately apply the Harmonic Oscillator Model of Aromaticity (HOMA)^55,56 as a structure-based aromaticity index under pressure, reference values for a C–C single bond and a C=C double bond were first determined. Ethane and ethene were selected as the smallest possible molecules for this purpose. Pressure calculations were performed using the GOSTSHYP method at the B3LYP-D3/cc-pVTZ level of theory, and the carbon–carbon bond lengths were extracted at each pressure. These reference values enable accurate evaluation of HOMA under compression by accounting for pressure-induced changes in the reference parameters.

To determine Nucleus-Independent Chemical Shifts (NICS) indices⁵⁷ based on pressure-optimized structures of benzene in both gas phase and crystalline environments, NMR calculations at the B3LYP-D3/cc-pVTZ level of theory were performed using Q-Chem. Ghost atoms were placed at the ring centers and 1 Å above them. Additionally, line scan analyses were conducted and fitted to models of infinitely-thin conducting lines of current.⁵⁸ Moreover, Gauge-Including Magnetically Induced Current (GIMIC) calculations were performed to directly evaluate the ring current at varying pressures.⁵⁹ The perturbed density matrices for this purpose were computed using Q-Chem 6.0 and Gaussian16.⁶⁰ To minimize numerical artifacts, the average over the six C–C bonds was taken.

For the determination of Ring Critical Point (RCP) properties, wavefunction files were generated for the optimized benzene structures in both environments using Q-Chem, followed by Quantum Theory of Atoms in Molecules (QTAIM)⁶¹ analyses with Multiwfn 3.8.⁶² Delocalization indices were computed with the AIMAll program⁶³ and multicenter indices were evaluated with a combination of the AIMAll, APOST-3D⁶⁴ and ESI-3D⁶⁵ programs.

Results

Structural indices

Since application of hydrostatic pressure to an isolated benzene molecule via X-HCFF and GOSTSHYP leads to isotropic compression, no changes in the symmetry of benzene can be observed. Hence, structural indices cannot be used feasibly to assess pressure-induced changes in aromaticity of benzene.

However, benzene crystallizes in an orthorhombic structure with the space group P_bca. Within the crystal, the molecules are arranged in a layered fashion and stabilized by van der Waals forces, quadrupole interactions and weak C–H⋯π hydrogen bonds between adjacent aromatic rings.^66,67 Due to the anisotropic environment within the crystal lattice, an individual benzene molecule experiences slight distortions compared to its isolated gas phase counterpart. Specifically, a subtle bond-length alternation arises, resulting in non-equivalent C–C bond lengths within the benzene ring. This deviation from the perfect hexagonal symmetry leads to a minor reduction in aromaticity, as the uniform delocalization of π-electrons is subtly perturbed. To further investigate pressure-induced changes in aromaticity within the crystal, various geometry-based indices were employed.

Bond-length alternation (BLA)

The Bond-Length Alternation (BLA) index is one of the most straightforward and commonly used structure-based indices. It quantifies the variation in bond lengths within a molecule serving as an indicator of electronic conjugation. BLA is typically understood as measuring the difference between the lengths of the longer single and shorter double C–C bonds in a π-conjugated system.^68,69 A lower BLA value corresponds to more uniform bond lengths, reflecting stronger delocalization, while a higher BLA indicates enhanced localization and reduced aromaticity. However, several studies have highlighted that π-electrons can exhibit a tendency to localize double bonds, counteracting the delocalizing influence of σ-electrons.^70,71 Accordingly, structural descriptors should be applied with caution when quantifying aromaticity.^12,72

In the case of an asymmetric benzene molecule within the crystal structure, BLA index can be calculated as the difference between the longest and shortest C–C bond lengths R within the aromatic ring.

BLA (def1) = R_max − R_min (2)

To account for the full atomic sequence of the benzene ring, BLA can alternatively be defined as the absolute difference between the average lengths of the even- and odd-numbered C–C bonds along the ring.^62,73 This approach provides more detailed insight into how uniformly bond-length alternation is distributed across the ring, especially in anisotropic environments such as the crystalline phase under pressure.

In both definitions, an increase in BLA signifies reduced structural aromaticity, as it reflects greater localization of the C–C bonds within the ring. For reference, benzene in its ideal gas phase form has a BLA of 0, reflecting its fully equalized C–C bond lengths and highly symmetric aromatic character. Regarding the first definition, the evolution of BLA under pressure up to 100 GPa is shown in Fig. 1.

Fig. 1 BLA index (def1, eqn (2)) for a single benzene molecule in the crystal lattice (Crystal Geo) from 0 to 100 GPa.

At ambient conditions, the BLA value is slightly above zero as expected due to the inherent anisotropy in the crystal environment. As pressure increases, a linear rise in BLA is observed between 5 and 25 GPa, followed by a gradual saturation at higher pressures. This trend reflects the increasing proximity of molecules under compression, which enhances intermolecular interactions and amplifies the anisotropic pressure acting on each benzene molecule. These forces lead to structural distortions and a gradual increase in bond-length alternation. The structural distortions appear to reach a threshold, beyond which further pressure has a diminishing effect. This behavior is consistent with the expectation that, at very high pressures, molecular packing becomes extremely dense, restricting additional structural rearrangements. Overall, the observed changes in BLA remain relatively small, suggesting that, despite pressure-induced structural distortions, benzene largely retains its aromatic character within this pressure range.

In contrast, BLA (def 2) does not reveal a clear trend under pressure, as the observed changes are extremely small, on the order of 10⁻⁶ Å (Fig. S4). This discrepancy highlights that different BLA definitions rely on distinct symmetry reductions and therefore yield different results. To complement these BLA measures, we introduce a self-defined parameter accounting for the average difference between all C–C bond lengths at each pressure. In the case of benzene, which contains six C–C bond lengths, all pairwise comparisons yield 15 possible differences, from which the mean value is subsequently determined. The resulting trend (Fig. S5) closely follows that of the BLA index (def 1), substantiating that bond-length alternation increases with pressure.

Harmonic oscillator model of aromaticity (HOMA)

Another widely used geometry-based descriptor of aromaticity is the HOMA index.^9,55 HOMA evaluates each of the n C–C bonds within an aromatic ring by comparing the length of each C–C bond i to a reference value Rref. This optimal bond length is selected such that the compression energy of the double bond and the expansion energy of the single bond are minimal according to the harmonic potential model.⁷⁴ The normalization factor α ensures that, for localized double and single bonds as defined by the harmonic potential model, HOMA equals to zero.^55,56

From this definition, a fully aromatic system in which all C–C bond lengths equal to R_ref exhibits a HOMA value of 1. In contrast, nonaromatic compounds typically display HOMA values close to 0, while antiaromatic monocyclic compounds often exhibit negative HOMA values. In this sense, HOMA measures the geometrical similarity to archetypal aromatic compounds⁷⁵ and fail in analyzing aromaticity changes along the Diels–Alder reaction.⁷⁶ A reparametrized version of HOMA has recently been introduced in which the empirical parameters have been derived from high-level quantum chemical calculations using both aromatic and antiaromatic reference compounds.⁷⁷ Nevertheless, the number of situations for which reference values have been tabulated remains limited, restricting the broader applicability of HOMA in inorganic compounds.

To determine appropriate values for R_ref and α under pressure, ethane and ethene were used as reference compounds. The bond lengths of their central single and double bonds were measured across a range of pressures. Detailed information on these calculations, along with tabulated values of R_ref and α and the corresponding bond lengths at different pressures, is provided in the SI (Table S2). Fig. 2 illustrates the pressure-dependent evolution of the HOMA index.

Fig. 2 HOMA index using ethane and ethene as reference systems for a single benzene molecule in the crystal lattice (Crystal Geo) from 0 to 100 GPa.

At zero pressure, the HOMA value is close to 1, indicating that the anisotropic environment of the crystal lattice has a negligible impact on the molecular geometry of individual benzene molecules. Aside from a slight initial increase, HOMA gradually decreases up to 55 GPa. This trend is consistent with the behavior observed in the BLA index, further reinforcing the conclusion that increasing pressure induces structural distortions that progressively weaken the π-electron delocalization within the benzene ring. However, the relative changes in HOMA are exceedingly small, suggesting that the overall variation in aromaticity remains minimal.

In addition, butadiene was tested as an alternative reference system. When these reference parameters are employed in the HOMA calculation, a different pressure-dependent curve trend is obtained as compared to ethane and ethene (Fig. S9). This highlights the importance of the reference parameters in HOMA-based analyses, especially when the observed changes are rather small, as in the case for benzene. In such instances, the choice of reference system can significantly influence the results and may lead to different or even contradictory trends.

In contrast to BLA, the evolution of the HOMA index can also be determined for the fully symmetric gas phase molecule with identical C–C bond lengths. For this purpose, pressure-optimized molecular geometries obtained with GOSTSHYP were employed. Fig. 3 shows the evolution of HOMA using ethane and ethene as reference systems.

Fig. 3 HOMA index using ethane and ethene as reference systems for a benzene molecule in the gas-phase optimized with GOSTSHYP (Gas – GOSTSHYP) from 0 to 100 GPa.

Similar to the crystalline molecule, the index exhibits an initial value close to unity followed by a minimal increase of up to 5 GPa, and subsequently decreases over the remaining pressure range. Although the variation is more pronounced than in the crystal phase, the overall change remains limited. When butadiene is used as the reference system for the gas-phase molecule, an almost linear decreasing trend is observed (Fig. S11). Employing ethane and ethene as reference systems yields more intuitive values, particularly at 0 GPa. This observation can likely be attributed to the partial delocalization of the four π-electrons in butadiene, which consequently affects the reference bond lengths used in the HOMA calculation.⁷⁸

Additionally, the evolution of the HOMA index was examined without pressure parameterization for both the crystal molecule (Fig. S10) and the isolated molecule (Fig. S12). In both cases, a pronounced decrease was observed even at low pressures, with this effect being more significant for the crystal molecule. At high pressures, the HOMA index decreased to ranges that no longer permit meaningful interpretation. This underscores the necessity of employing appropriate reference values under pressure conditions.

Energetic indices

The concept of describing aromaticity through an energetic measure originates from Pauling's valence bond theory.⁷⁹ One of the hallmark features of aromatic systems is the additional stabilization conferred by the delocalization of π-electrons within a cyclic conjugated system. Numerous methods have been proposed to quantify the so-called Aromatic Stabilization Energy (ASE).⁸⁰ The main challenge in the quantification of ASE lies in the selection of an appropriate nonaromatic reference system for the molecule of interest. Ideally, this reference should be carefully selected to eliminate contributions from other effects, such as strain, hyperconjugation and hybridization changes, thereby isolating the elusive contribution from aromaticity.⁸ To distinguish between different types of reactions used to evaluate ASE, a classification into isogyric, isodesmic, hypohomodesmotic, homodesmotic and hyperhomodesmotic reactions can be made based on the conservation of bond and atom types.⁸¹ In this study, two reaction schemes are employed, involving a hypohomodesmotic reaction and a homodesmotic reaction (Fig. 4).⁸² The homodesmotic reaction offers the advantage of eliminating spurious contributions from hyperconjugation and hybridization effects, but it might introduce artifacts due to the involvement of multiple molecular species.^80,83 At zero pressure, the ASE derived from the homodesmotic reaction is estimated to be −34.1 kcal mol⁻¹ at the B3LYP-D3/cc-pVTZ level of theory, which is in excellent agreement with ASE values previously reported in the literature for benzene.^8,84

Fig. 4 Evolution of the aromatic stabilization energy of benzene as a function of pressure, derived from two distinct reaction schemes: (a) a hypohomodesmotic reaction and (b) a homodesmotic reaction.

The evolution of the ASE as a function of pressure for the two distinct reaction schemes is displayed in Fig. 4. In these reactions, all species are compressed at the same level of theory. It is evident that the trends in the stabilization energies are highly dependent on the considered reaction. The homodesmotic reaction on the right is only weakly impacted by geometric distortion (X-HCFF), but suggests an increase in aromatic stabilization of about 5 kcal mol⁻¹ at 100 GPa due to electron confinement (GOSTSHYP). In contrast, the hypohomodesmotic reaction on the left exhibits the opposite behavior: ASE increases upon geometric contraction, while it decreases upon electronic confinement. It should be noted that the absolute value of the volume reduction/expansion in common chemical reactions can easily reach ≈10 cm³ mol⁻¹.⁸⁵ This will result in enthalpy corrections that can exceed the observed changes in aromatic stabilization energies by two orders of magnitude. Thus, while the observed energetic variations are significant relative to gas-phase values, they remain consistent with a limited variation in aromaticity under high-pressure conditions.

A reason behind the opposing trends for the two reaction schemes can be better understood by considering the first-order change in the reaction energy ΔE_R in eqn (5). Note that the f and i subscript denote the collection of products (final) and reagents (initial), respectively.

here, the identities for pressure and compressibility are used. The first-order change in ΔE_R clearly depends on the compressibility (which is in turn influenced by the force constants of highly delocalized skeletal vibrations⁸⁶) as well as the volume, which depends not only on the molecular geometry but also on the connectivity. Neither of these effects are accounted for by the reaction classes described above, beyond the local level. As a result, the effect of pressure on aromaticity, as mediated through the compressibility, is difficult to isolate robustly. A meaningful analysis would require a process that does not significantly change the system's geometry, while a pragmatic approach could match geometrical changes with the application of interest. However, in practice, the enthalpic contribution from volume changes dominates, explaining the prevalence of reaction and activation volumes in high-pressure chemistry literature.⁸⁵

Electronic indices

Ring critical point (RCP) properties. The delocalization of electrons within a ring can be assessed through Atoms In Molecules (AIM) parameters.^61,87 In this framework, aromaticity is generally evaluated at the Ring Critical Point (RCP) based on several quantities being the electron density, its Laplacian, and the total energy density H along with its components, namely the kinetic energy density G and potential energy density V:^88,89

H = G + V (6)

Typically, regions with greater electron delocalization are characterized by high electron densities, Laplacians with a positive sign, and relatively large positive values of G. In contrast, regions where electrons are more localized exhibit low electron densities, negative Laplacians, and relatively large absolute values of V (negative). Furthermore, the total electronic energy H is positive, as G dominates over V. In phenyl rings, an increase in both electron density and its Laplacian is accompanied by an increase in total energy density H and the absolute values of its components, G and V, indicating an increase in aromaticity.⁸⁹

Fig. 5(a)–(e) illustrates the pressure-dependent evolution of AIM-based parameters for benzene, including the electron density, the Laplacian, as well as the energetic descriptors (H, G, and V) evaluated at the RCP. Fig. 5(f) shows the average distance of each carbon atom to the geometric center of the ring for the gas-phase molecule under compression using X-HCFF and GOSTSHYP as well as for a benzene molecule in the crystalline phase.

Fig. 5 Pressure-dependent evolution of AIM descriptors at the RCP for gas phase and crystal benzene using different compression models: (a) electron density, (b) Laplacian of electron density, (c) total electron energy density H, (d) kinetic electron energy density G and (e) potential electron energy density V at the RCP. In (f), the mean distance of each carbon atom to the geometric ring center for a benzene molecule in the gas phase and in the crystal is shown. The abbreviations stand for the following: gas – X-HCFF: gas-phase molecule optimized using X-HCFF, gas – GOSTSHYP: gas-phase molecule optimized with GOSTSHYP, gas – GOSTSHYP Geo: single-point calculation without pressure of GOSTSHYP-optimized gas-phase geometry at each pressure, Crystal Geo: single-point calculation on a molecule extracted from the crystal molecule at each pressure, Crystal Geo – GOSTSHYP: single-point GOSTSHYP calculation on a crystal-extracted molecule at each pressure.

All descriptors for both the gas phase and crystal molecules begin at nearly identical values, indicating that the slight asymmetry in the crystal benzene structure has only a minor influence on the electronic properties at the RCP under ambient conditions. Upon compression of the ring, both the electron density and its Laplacian at the RCP increase, accompanied by rising absolute values of G and V. In this regime, G dominates over V, which generally yields positive values for H. However, except for GOSTSHYP, the change in absolute potential energy under pressure is more pronounced than that in kinetic energy, resulting in decreasing trends for H. This behavior leads to seemingly contradictory interpretations of the evolution of aromaticity under pressure.

The pressure-dependent aromaticity trends also vary with the compression method. Benzene experiences the highest compression with X-HCFF, followed by the crystal molecule, while GOSTSHYP produces the weakest compression (Fig. 5f). Consequently, the X-HCFF compressed molecule shows markedly steeper variations in all indices compared to both the crystal molecule and the GOSTSHYP results. Interestingly, the additional consideration of electron compression using GOSTSHYP has only a very modest impact on the electron density and its Laplacian at the RCP. In contrast, a stronger dependence is observed for the total electron energy density. While H decreases significantly under X-HCFF, it remains nearly constant across the entire pressure range with GOSTSHYP. Regarding the crystal molecule, a discontinuous trend is observed without compression of the electron density.

Overall, most of the indices, except for the total energy density H, point to a modest increase in aromaticity for both the gas phase and crystal molecules. Caution is advised, however, as compression of the ring inherently increases the electron density at the RCP, its Laplacian, and energy components G and V. These effects are not unique to benzene and can likewise be observed in nonaromatic systems such as cyclohexane (Fig. S38).

Delocalization indices. Another important class of electronic aromaticity descriptors involve the so-called Electron Sharing Indices (ESIs), which quantify the extent of electron delocalization between atoms or molecular fragments.^10,90 These indices are typically derived from the exchange-correlation density, defined as the difference between the real pair density and the fictitious pair density of non-interacting one-electron densities as given by eqn (7).^91,92 The pair density can be evaluated at any desired level of theory; however, corrections beyond a single-determinant approximation are generally small.⁹³

When electrons do not interact, the probability of finding an electron at is not influenced by the presence of an electron at . In this case, the exchange-correlation density vanishes and any nonzero value thus reflects the degree of electron–electron interaction. When integrated over two spatial regions, ρ_XC provides a measure of the extent to which electrons are delocalized across these regions. For this reason, ρ_XC has been frequently used in the definition of ESIs. Among these, the Delocalization Index (DI) δ(A, B) is used in the present work.^94–96 This index is obtained by integrating ρ_XC over the atomic basins of atoms A and B (eqn (8)), most commonly defined within the QTAIM framework from the topology of the electron density. Nevertheless, alternative atomic partitions may also be employed, such as Hilbert-space partitions, which provide reliable and computationally efficient alternatives to real-space partitioning.^93,97

Since δ(A, B) quantifies the number of electron pairs shared between two atomic regions, several aromaticity indices have been developed based on this concept.^10,90 The Para-Delocalization Index (PDI) is the average of the DIs for para-related atoms in six-membered rings, thereby providing a measure of the electron delocalization across the ring. Poater et al. demonstrated that this index correlates well with other established measures of aromaticity in polycyclic aromatic hydrocarbons.^98,99 A larger PDI value reflects greater aromaticity, but is restricted to six-membered rings and overestimate the delocalization in heterocycles.¹⁶

Another electronic index is the Aromatic Fluctuation Index (FLU), which measures the degree of delocalization between adjacent atoms along the ring and its deviation from an ideal aromatic reference value (δ_ref(A_i, A_i+1)).¹⁰⁰ Unlike HOMA (vide supra), this reference value is more difficult to be determined under pressure, as it must be derived from a suitable aromatic molecule. For C–C bonds, the aromatic reference is benzene. The expression for the FLU is given in eqn (9), where the fraction of atomic localization indices reflects the electron sharing similarity between δ(A_i) and δ(A_i+1) adjacent atoms, α is either 1 or -1 to ensure the fraction is larger than 1 and n is the number of adjacent atoms. FLU approaches zero for aromatic species and can be applied to rings of any size. However, since it measures the similarity with respect to a predefined aromatic molecule, FLU is not well suited for describing reactivity.⁹⁰

From the delocalization indices in Fig. 6, it is evident that geometric contraction and electronic confinement exert completely different effects. The contraction of the geometry results in an increased orbital overlap (or reduced interelectronic distances), thereby increasing the delocalization index compared to zero-pressure conditions. This effect diminishes with increasing C–C distance, i.e. from ortho-to meta-to para-related carbon atoms, and can even reverse under high pressure when significant geometric distortions occur. By contrast, incorporating electronic confinement through the GOSTSHYP model produces the opposite effect, indicative of electron localization at all pressures. In some cases, this electron localization can outweigh the effect of geometric contraction, as observed for the full GOSTSHYP calculation (Gas – GOSTSHYP) in the ortho- and para-positions, though this is not always the case. Additional factors, such as bond-length alternation induced by the crystal environment or adopting the Bondi set of atomic radii for the cavity construction, appear to have little influence on the delocalization indices (see Fig. S15 and S16), but summing the hydrogen atoms into the carbon basins decreases the effect of confinement somewhat, further indicating that confinement influences the structure of the π system.

Fig. 6 Change in the delocalization index at 50 and 100 GPa compared to 0 GPa for the ortho-, meta- and para-position using: (a) the X-HCFF method; (b) the GOSTSHYP method; (c) the GOSTSHYP geometries without electronic confinement and (d) the crystal-phase geometry combined with the GOSTSHYP method for electronic confinement.

These trends in the DIs translate into predictable changes in the PDI and FLU electronic indices. The FLU decreases under geometric contraction and increases under electron confinement (provided the reference parameters remain unaltered), whereas the PDI increases under geometry contraction and decreases with electron confinement. The increase of the PDI under geometry contraction was also observed by Feixas et al. in the context of aromaticity for distorted benzene.¹⁰¹ An interesting comparison is obtained by examining the PDI of the nonaromatic cyclohexane under pressure (SI). Its PDI also increases with pressure, as expected from the reduced interelectronic distances and enhanced electronic interactions, and this increase is even more pronounced than in benzene using the X-HCFF model. Taken together with the decrease in DIs observed under the full GOSTSHYP method, this behavior points to a diminishing electronic aspect of aromaticity under pressure. It should be noted, however, that the overall magnitude of the pressure-induced changes in electronic aromaticity are far less significant than those observed between different aromatic molecules.

Multicenter indices. While the delocalization indices quantify the degree of electron delocalization starting from pairs of atoms, another class of electronic descriptors generalizes the concept of bond orders from pairs of atoms to any number of atoms M.¹⁰² This lead to the definition of the multicenter bond order in terms of the density matrix P in eqn (10) and the working eqn (11), where S_ij(A) are the overlap integrals of molecular orbitals i and j in the region associated to atom A. These have been proposed as a measure of aromaticity when applied to the sequence of atoms A₁ to A_M forming a conjugated ring system.¹⁰³

This quantity can also be interpreted as the charge susceptible to delocalization along the ring, reminiscent of the magnetic criterion of aromaticity,¹⁰⁴ or as the possibility to form sequences of resonating orbital overlaps between neighboring atoms. A generalization of the I_ring is possible by removing the constraint of a circular atom sequence from A to M and instead summing over all possible path permutations, resulting in the multicenter index (MCI) introduced by Bultinck et al.¹⁰⁵ This extension enables the incorporation of cross-ring delocalization in the description of aromaticity. MCI emerges as one of the most reliable indices in several benchmark tests designed to evaluate the performance of aromaticity descriptors in both organic and inorganic species.^16,19 Both electronic indices adopt larger positive values in aromatic species, while nonaromatic molecules yield values close to zero. Antiaromatic systems exhibit intermediate I_ring and MCI values, making their identification more challenging. Additional variants have been developed, such as AV1245,¹⁰⁶ particularly suited for larger rings by mitigating numerical error propagation, and AV_min,¹⁰⁷ capable of detecting localized disruptions in the π-conjugated system.

For benzene under pressure, all four multicenter indices considered (i.e. I_ring, MCI, AV1245 and AV_min) consistently indicate a decrease in aromaticity with increasing pressure (Fig. 7a and b for I_ring and MCI, respectively). This is consistent with the reported decrease in multicenter indices observed during the bond-length compression of benzene.¹⁰¹ This effect appears to be dominated by geometric changes, with only a minor contribution from the electronic confinement. Overall, the decrease in the aromaticity indices is modest and the different partitions agree reasonably well (Fig. 7c and d), apart from the Hirshfeld partition, which suggests that electronic confinement markedly increases the aromaticity under pressure (Fig. S19). Using van der Waals radii instead of isosurface-fitted radii in the pressure models does not lead to significant changes (Fig. S19). In cyclohexane, the effects of pressure on multicenter descriptors are even smaller. Taken together, these results confirm that according to electronic criteria, aromaticity decreases under pressure.

Fig. 7 (a) I_ring index as a function of pressure for benzene using the QTAIM partition. MCI index as a function of pressure for benzene with the QTAIM partition (b), the topological fuzzy Voronoi cell partition (c) and the Hirshfeld partition (d).

Linear response function (LRF). The (frequency-independent) linear response function is a central quantity in conceptual density functional theory,^108,109 which has been previously applied to quantify electron delocalization in general^110,111 and aromaticity in particular.^112–114 The conceptual DFT framework is a perturbational expansion of the electronic energy of a system in terms of its fundamental variables: the number of electrons N and the external potential . Its aim is to establish rigorous mathematical definitions for often vaguely defined chemical concepts. Within this framework, the linear response function is the second-order functional derivative of the electronic energy with respect to the external potential and can be computed at several levels of approximation.¹¹⁵ The expression for the simplest case, i.e. the Independent Particle Approximation (IPA), is given in eqn (12).here, are the molecular orbitals and ε their orbital energies, where a denotes virtual and i occupied orbitals. Sablon et al. showed that the atom-condensed LRF kernel could discern between aromatic and nonaromatic systems and they proposed the Para-Linear Response (PLR) as a new quantitative aromaticity index, inspired by the PDI (vide supra).¹¹² The PLR is defined as the average value of the atom-condensed LRF for para-related carbons in six-membered rings and was proven to correlate well with the PDI. Subsequent studies further established the LRF as a suitable indicator for σ-aromaticity, metalloaromaticity, and even non-equilibrium states.^112–114 The behavior of the LRF can intuitively be understood by considering that a perturbation at the ipso carbon leads to a more delocalized response, and consequently stronger interaction at other atoms within the ring, in aromatic systems. In this work, the LRF kernel is condensed to atoms by integration over the Hirshfeld atomic basins, and its value for different carbon–carbon positions along the aromatic benzene ring is employed as a measure of aromaticity.

In general, the linear response decreases at all positions with increasing pressure, consistent with a reduction in aromaticity. This trend is mainly driven by geometric distortions of the benzene ring and can be rationalized by the greater spread of the occupied and virtual orbital energies at shorter bond lengths. Confinement plays only a minor role, as evidenced by the similarity between the GOSTSHYP and GOSTSHYP Geo methods (Fig. 8a and b). Electronic confinement compresses the virtual orbitals, increasing their overlap with the occupied orbitals and exerting a slight positive effect. Following the relative magnitude of geometric distortions induced by the different models, the X-HCFF model (Fig. 8c) predicts a larger decrease in aromaticity than GOSTSHYP, while the crystal structure takes an intermediate position. No particular anisotropy effects are observed (Fig. S13). At 100 GPa, the PLR indicates a reduction in aromaticity of 10–20%. The LDA approximation (Fig. 8d) provides the same picture as the IPA approximation with similar relative changes, but at substantially higher computational cost (Fig. S14).

Fig. 8 Matrix elements of the linear response kernel for differently related carbon atoms along the benzene ring at varying pressures: (a) IPA approximation with the GOSTSHYP model; (b) IPA approximation with the GOSTSHYP Geo model; (c) IPA approximation with the X-HCFF model and (d) LDA approximation with the GOSTSHYP model.

Since the Hirshfeld partitioning is employed to condense the LRF, it is feasible to construct the partition function at high pressures from pressurized atomic densities. This leads to almost no changes in the PLR with pressure, but summing the H-atom regions into the C regions restores the trend without pressurized atomic densities (Fig. S13). Otherwise using CH regions has no effect on the observed effects of pressure.

Magnetic indices

Nucleus-independent chemical shift (NICS). NICS is determined by evaluating the magnetic shielding at any given point in space, thereby probing the strength of the ring current effects independently of the atomic nuclei. It is defined as the negative of the total magnetic shielding. Negative NICS values generally indicate aromaticity arising from diatropic ring currents, while positive values suggests antiaromaticity associated with paratropic ring currents.⁵⁷ The isotropic NICS(0) index is measured at the center of the ring in the molecular plane and thus includes contributions not only from the π-electrons, but also from unpaired electrons, core electrons, and the local distribution of σ-bond electrons. To minimize these additional effects, the isotropic NICS(1) value was introduced, measured 1 Å above the molecular plane.¹¹⁶ While the total magnetic shielding is isotropic, the shielding tensor also contains contributions from local currents induced by fields that are not perpendicular to the ring. These contributions are not representative for the aromatic ring current and are thus undesirable. A more precise description of electron delocalization can be obtained from the out-of-plane tensor components NICS(0)_zz and NICS(1)_zz.^117,118

Fig. 9 shows the evolution of the different NICS-based descriptors under pressure. The limited impact of crystal anisotropy on magnetic aromaticity under ambient conditions is evident, as all indices begin at nearly identical values. For NICS(0)_zz, NICS(1), and NICS(1)_zz, a consistent downward trend is observed, irrespective of the molecular environment or the compression method, indicating an increase in aromaticity. This agrees with the observed behavior of the magnetic descriptors computed along the symmetric ring stretching normal mode.^101,119 Conversely, NICS(0) displays an upward trend when electron density compression is not included. These upward deviations arise from in-plane magnetic responses, as evident from their absence in the evolution of the out-of-plane component NICS(0)_zz. Within the molecular plane, electron density compression via GOSTSHYP strongly affects the behavior of NICS(0) and NICS(0)_zz, whereas above the molecular plane its impact is notably lower as seen from the NICS(1) and NICS(1)_zz curves. In the absence of electron density compression, the most pronounced trends are obtained with X-HCFF, consistent with this method inducing the largest molecular compression (Fig. 5f).

Fig. 9 Evolution of (a) NICS(0), (b) NICS(1), (c) NICS(0)_zz and (d) NICS(1)_zz for benzene from 0 to 100 GPa.

In general, confinement of benzene under pressure leads to more negative NICS-based indices, consistent with an increase in aromaticity. Nevertheless, even for NICS indices known to overestimate aromaticity,¹²⁰ the magnitude of these changes remains relatively small.

Infinitely-thin conducting line of current (ICLOC) model. Despite their usefulness and wide application, NICS-based indices remain an indirect probe of the current flow associated with aromaticity.³ Ring current models, introduced in the 1950s, provide a more direct way to describe current susceptibility through the ring current parameter that is more closely associated with aromaticity.^121–123 These models simplify the quantum-mechanical current density into a few meaningful parameters, by representing it as classical electromagnetic currents. The simplest is the Infinitely-thin Conducting Line Of Current (ICLOC) model,¹²⁴ in which the total current induced by a magnetic field is carried through a single loop with radius R. More elaborated models combine loops positioned in, above or below the molecular plane and carrying either a diatropic or a paratropic current.^58,124–126 Alternatively, the loop can be assigned a finite thickness, resulting in a toroidal model (TCLOC), or a polygonal shape.^127,128 The parameters of these models are typically determined by fitting the magnetic shielding they produce to the NICS(0)_zz values along the axis perpendicular to the molecular plane through the ring center.

As the size of the benzene ring decreases under pressure, the radius of the ring current loop is expected to decrease accordingly, leading to an increase in the magnetic shielding that is unrelated to the ring current typically associated with aromaticity.^129,130 To investigate this effect, the quantum-mechanically calculated NICS values were fitted to their classical counterparts from different ICLOC models (Fig. 10). The considered models include a single loop of current (ICLOC1), a single toroidal loop of current (TCLOC1), a double loop above and below the plane (ICLOC2Z), a double diatropic and paratropic loop in the plane (ICLOC2XY), a combination of these two double-loop models (ICLOC4) and a hexagonal variant of the latter model (ICLOC4hex). Since the NICS values at different pressures along the line scans differ only slightly, high-quality fits were essential to extract meaningful parameters. To achieve this, several fitting procedures were employed, including a standard least-square fit (LSQ), a least-squared fit up to 5 Å (LSQ5A), a weighted least-squares fit with the square of the NICS as the weighting factor (WSQY2) and a weighted least-squares fit using the variance of the NICS between different pressures as the weighting factor (WSQVAR).

Fig. 10 Schematic illustration of four ICLOC models corresponding to a single infinitely thin loop (ICLOC1), two elevated loops above and below the plane (ICLOC2Z), a single loop with a finite thickness (TCLOC1), two loops in the plane carrying a paratropic and diatropic current (ICLOC2XY) and two loops above and below the plane and two loops in the plane carrying a paratropic and diatropic ring current (ICLOC4).

Upon inspection of the fitted models, the ICLOC1 model was immediately discarded as it fails to reproduce the reduced response within the plane of the benzene ring. We further imposed the requirement of obtaining a successful fit for every pressure, which often proved problematic for the seven-parameter models (ICLOC4 and ICLOC4hex). In addition, visual inspection of the region where pressure effects are most pronounced (i.e. ≈1 Å) eliminated models that were unable to capture the distinctive features in this region. Full details on all models and selected fitting parameters are provided in the SI (S22–S33).

The trends in the parameters of the remaining models, after excluding the least suitable ICLOC variants, showed rather good agreement. The TCLOC1 model proved to be particularly suitable, while the ICLOC4 and ICLOC2XY models also performed adequately in selected cases, occasionally requiring the WSQY2 fitting procedure.

For all pressure models, the ring current parameter (I) increases with pressure (Fig. 11), consistent with the rise in magnetic aromaticity predicted by NICS-based indices. Among these indices, the NICS(1)_zz best reproduces the overall ring-current trends; however, at 100 GPa, the ring current increases by only 2–4%, roughly half the increase suggested by the NICS(1)_zz. The remaining increase arises from the shrinking radius of the current loop, which decreases proportionally to the interatomic distances and mirrors the bond compression observed for the different pressure models. When included in the model, the Z parameter controlling the elevation of the current loops above and below the plane follows the same trend. Both the enhanced ring current and the reduced current loop radius can be attributed to the geometric distortion and the electronic compression. For crystal geometries at the highest pressures, however, the ring current decreases again, indicating a loss of aromaticity that is not fully captured by the NICS indices. In the TCLOC1 model, the additional A parameter controls the width of the current band, which broadens under geometric contraction under pressure but narrows when electron density compression is also included, although this effect is minor. The same analysis with the TCLOC1 model and LSQ fitting procedure using the isosurface-fitted radii, provided excellent agreement with these trends.

Fig. 11 (a) Ring current and (b) radius parameter of the toroidal line of current model as a function of pressure with various models for benzene.

Gauge including magnetically indiced current (GIMIC). Analogous to ICLOC models, the electronic flux can directly probe the ring current related to aromaticity and is thus able to show the most subtle effects.^131,132 Because of the origin dependence of the magnetic potential , the vector field also depends on the chosen origin in a finite basis. The use of Gauge-Including Atomic Orbitals (GIAOs) minimizes this dependence, rendering the flux gauge-independent, though not fully invariant,¹³³ as implemented in the GIMIC program.^134,135

The ring current strength can be quantified by integrating the flux over the plane bisecting a C–C bond, defined by symmetry as the plane passing through the geometric center of the benzene ring and the bond midpoint.^135,136 Since the integrated current density profiles obtained from Q-Chem calculations exhibited nonphysical behaviour, only the Gaussian16 results are reported (Table 1). Since these Gaussian calculations do not explicitly include confinement effects, the ICLOC analysis provided above remains valuable as a complementary descriptor of magnetic aromaticity.

Table 1 GIMIC ring current strengths (nA/T) for pressure-distorted benzene geometries, calculated from perturbed density matrices obtained with Gaussian16

	Single – GOSTSHYP Geo		Crystal Geo
	Dia	Para	Dia	Para
0 GPa	17.0	−5.0	16.7	−4.9
20 GPa	17.2	−5.0	17.0	−5.0
40 GPa	17.4	−5.0	17.1	−5.0
60 GPa	17.6	−5.0	17.2	−5.0
80 GPa	17.7	−5.0	17.2	−5.0
100 GPa	17.8	−5.0	17.2	−5.1

---

These results indicate that the ring current strengthens under pressure as the diatropic component increases. This finding is consistent with previous studies based on the magnetic criterion and points to an enhancement of aromaticity associated with the geometric distortion induced by pressure.

Comparison of aromaticity indices and recommendations. To compare the various aromaticity descriptors, the relative trends of all indices with respect to their maximum change are shown at pressure intervals of 10 GPa (Fig. 12). For clarity, these trends are presented exclusively for the GOSTSHYP model applied on a single molecule (except for the BLA and HOMA which consider the crystal structure), as this model provides a quantum-mechanical description of both geometric and electronic effects induced by pressure. The corresponding trends of all indices for both the gas phase and the crystal molecule across different compression models are provided in Fig. S35 and S36.

Fig. 12 Normalized evolution in aromaticity under pressure at 10 GPa intervals. For each index, normalization is performed with respect to the maximum observed change in aromaticity over the whole pressure range, using the value at 0 GPa as the reference (yellow). An increase in aromaticity is represented by red coloring, whereas a decrease is indicated by green coloring. Descriptors where decreasing values indicate increasing aromaticity are shown above the bold dividing line; while those where increasing values correspond to greater aromaticity are shown below. For the energetic, magnetic and electronic descriptors, only the gas phase GOSTSHYP results are displayed. White boxes are inserted where no data are available (see full computational details in the SI).

Overall, only minor changes are observed across all indices up to 100 GPa, regardless of the underlying criterion, suggesting that electron delocalization is largely preserved even at high pressure. The most profound changes are found for the energetic and magnetic indices, while the electronic and structural descriptors display only subtle variations. Notably, a dichotomy emerges under pressure: magnetic descriptors consistently indicate an increase in aromaticity, while the intrinsic descriptors (i.e. electronic, structural and energetic) instead suggest a decrease. It should be noted, however, that the energetic descriptors are not internally consistent. Taken together, this divergence underscores the need for a diverse set of indices, since no single descriptor fully captures the multifaceted nature of aromaticity.

Structure-based aromaticity indices generally predict a slight decrease in aromaticity under pressure for both isolated and crystal benzene molecules. The BLA (def 1) accounts for only two bond paths within the ring, offering limited insight into structural changes, while the BLA (def 2) is restricted by symmetry. To address these limitations, an alternative definition was proposed based on the mean difference across all six C–C bond lengths, which reproduced the same trend as BLA (def 1). The HOMA index was reparameterized for non-ambient conditions to obtain pressure-dependent reference parameters. Ethane/ethene and butadiene were used as reference systems. The evolution of HOMA under pressure strongly depends on the chosen reference system, with ethane/ethene being recommended over butadiene.

Energetic indices are shown to be highly dependent on the considered reaction. Their application to evaluate resonance energies and aromatic stabilization energies under pressure is ambiguous and warrants further exploration. Additional factors, such as compressibility and reaction volumes, further complicate the isolation of the elusive contribution from aromaticity.

The RCP properties as electronic descriptors predominantly suggest an increase in aromaticity under pressure. However, similar trends in the electron density, its Laplacian, and the kinetic and potential energy densities components are also observed in cyclohexane, a nonaromatic molecule, indicating that these changes are a direct consequence of electronic structure compression rather than aromaticity. Only the total energy density H exhibits distinct behavior for the two molecules under GOSTSHYP compression, decreasing in cyclohexane while remaining nearly constant in benzene, thus supporting only minor changes in aromaticity under pressure. Nonetheless, RCP-based descriptors alone are insufficient to provide a definitive conclusion regarding the pressure-induced evolution of benzene's aromaticity; a limitation further compounded by their previously demonstrated unreliability in heterocyclic compounds.¹³⁷ Delocalization indices suffer from comparable limitations, though condensed values such as PDI offer clearer trends. By contrast, multicenter indices display consistent behaviour across all definitions, with the popular QTAIM partition remaining the recommended one.

NICS(1)_zz, a generally reliable magnetic index for monocyclic systems, shows a clear increase in aromaticity under pressure. Similar trends are found for NICS(0), NICS(0)_zz and NICS(1) when employing GOSTSHYP (Fig. 9), but the strong influence of pressure on the in-plane tensor components discourages the use of isotropic NICS(0) and NICS(1). Importantly, magnetic indices probe the response to an external magnetic field and therefore do not directly describe intrinsic aromatic properties under pressure.^13,138 In cases of rings with different compressibilities, an ICLOC approach is suggested to compensate for the decreasing ring size. Otherwise, the qualitative effects are well captured by the NICS(1)_zz.

In general, the different compression models yield broadly consistent trends, with variations largely attributable to differences in the predicted geometry compression and/or inclusion of electronic confinement. To quantify aromaticity under pressure, a combination of descriptors based on different criteria remains necessary and should include magnetic descriptors. While some indices, such as structural, multicenter and LRF, are relatively insensitive to additional electronic compression, others, including ASE, DI and NICS, are affected and must account for this confinement effect. Ultimately, a pressure model that integrates both geometric and electronic effects is essential for quantifying aromaticity, a recommendation likewise emphasized in the modeling of high pressure chemistry.³³

Conclusion and outlook

This work addresses the fundamental question of how aromaticity can be reliably evaluated under pressure by examining the evolution of widely used descriptors in benzene as the archetypal test case. A broad set of structural, energetic, electronic, and magnetic indices were critically assessed, revealing that the overall changes in aromaticity up to 100 GPa remain modest. While structural and electronic indices consistently indicate a slight decrease of aromaticity under pressure, all magnetic descriptors point to an enhanced ring current, underscoring the need for a multidimensional framework that integrates complementary descriptors.

Indices relying on reference parameters derived at ambient conditions, such as HOMA and FLU, may require reparameterization to yield meaningful results. Energetic descriptors are strongly dependent on the chosen reaction and reference molecules, and ASEs are perturbed by additional effects such as compressibility and reaction volumes. Isotropic NICS-based indices contain spurious in-plane tensor contributions; therefore, the use of the NICS(1)_zz, ICLOC models and ring current analysis is recommended to assess ring current strength under pressure. Multicenter indices do not rely on external parameters and provide a consistent description of pressure-induced aromaticity changes. Since the different indices differ in their sensitivity to geometric distortion and electron confinement, pressure models that simultaneously account for both effects at the molecular wavefunction level are essential for a reliable quantum-mechanical description of aromaticity under pressure.

Due to the rigidity of benzene, its aromaticity is largely preserved under compression, as evidenced by the limited extent to which the different aromaticity indices change under pressure. Building on this observation, we propose that the response of a molecular system to pressure may itself serve as an informative indicator of aromaticity. Highly aromatic systems, such as benzene, are expected to exhibit only minor pressure-induced changes in geometry, electron density, and aromaticity descriptors, whereas antiaromatic systems may display more pronounced pressure-induced variations in these indices. To assess this hypothesis, future work will extend the present analysis to larger polycyclic systems, such as acenes and helicenes, as well as to antiaromatic molecules (for which some preliminary data is included in the SI), where more substantial structural and electronic distortions under pressure are anticipated. Such studies will also help identify the descriptors that most reliably capture pressure-induced changes in aromaticity.

Author contributions

J. E. – conceptualization, methodology, software, validation, analysis, writing – original draft, review & editing; A. D. – conceptualization, methodology, software, validation, analysis, writing – original draft, review & editing; A. P. – methodology, software, validation, analysis; F. Z. – software; T. N. – conceptualization, methodology, resources, analysis, writing – original draft, review & editing, supervision; F. D. P. – conceptualization, resources, writing – review & editing, supervision; M. A. – conceptualization, methodology, resources, analysis, supervision, writing – original draft, review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

Most of the data supporting this article have been included as part of the supplementary information (SI). Supplementary information: full details of the pressure protocols (X-HCFF, GOSTSHYP, and solid-state benchmarks) and the aromaticity descriptors (structural, electronic, magnetic, and energetic), including definitions, reparametrization, and validation across levels of theory and partitioning schemes, together with preliminary results for antiaromatic charged benzenes. See DOI: https://doi.org/10.1039/d5sc07920a.

Acknowledgements

F. D. P. and M. A. thank the VUB for the Strategic Research Program awarded to the ALGC research group. J. E. acknowledges the FWO for their financial support (1148522N). He also thanks Dr Bin Wang for his assistance in computing the linear response kernel. A. D. acknowledges funding by the Fonds der Chemischen Industrie. The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation – Flanders (FWO) and the Flemish Government.

Notes and references

1. G. Merino, M. Solà, I. Fernández, C. Foroutan-Nejad, P. Lazzeretti, G. Frenking, H. L. Anderson, D. Sundholm, F. P. Cossío, M. A. Petrukhina, J. Wu, J. I. Wu and A. Restrepo, Chem. Sci., 2023, 14, 5569–5576.
2. M. Solà, A. I. Boldyrev, M. K. Cyrański, T. M. Krygowski and G. Merino, Aromaticity and Antiaromaticity: Concepts and Applications, Wiley, 2022, pp. 1–294.
3. Aromaticity: Modern Computational Methods and Applications, ed. I. Fernandez, Elsevier, 2021, pp. 1–499.
4. H. Szatylowicz, P. A. Wieczorkiewicz and T. M. Krygowski, Sci, 2022, 4, 24.
5. H. Ottosson, Chem. Sci., 2023, 14, 5542–5544.
6. M. D. Peeks, Adv. Phys. Org. Chem., 2024, 58, 1–37.
7. G. Frenking and A. Krapp, J. Comput. Chem., 2006, 28, 15–24.
8. M. Alonso and I. Fernández, in Aromaticity, Elsevier, 2021, pp. 195–235.
9. T. M. Krygowski, H. Szatylowicz, O. A. Stasyuk, J. Dominikowska and M. Palusiak, Chem. Rev., 2014, 114, 6383–6422.
10. F. Feixas, E. Matito, J. Poater and M. Solà, Chem. Soc. Rev., 2015, 44, 6434–6451.
11. R. Gershoni-Poranne and A. Stanger, Chem. Soc. Rev., 2015, 44, 6597–6615.
12. F. De Vleeschouwer, E. Desmedt and M. Alonso, Chem. Meth., 2025, 5, e202500064.
13. J. Yan, T. Slanina, J. Bergman and H. Ottosson, Chem. –Eur. J., 2023, 29, e202203748.
14. M. Baranac-Stojanović, Chemistry, 2025, 7, 127.
15. M. Solà, Front. Chem., 2017, 5, 22.
16. F. Feixas, E. Matito, J. Poater and M. Solà, J. Comput. Chem., 2008, 29, 1543–1554.
17. C. Silva López and O. Nieto Faza, in Aromaticity, Elsevier, 2021, pp. 41–71.
18. M. K. Cyrański, T. M. Krygowski, A. R. Katritzky and P. v. R. Schleyer, J. Org. Chem., 2002, 67, 1333–1338.
19. M. Solà, F. Feixas, J. O. C. Jiménez-Halla, E. Matito and J. Poater, Symmetry, 2010, 2, 1156–1179.
20. R. Báez-Grez and R. Pino-Rios, ACS Omega, 2022, 7, 21939–21945.
21. T. Stauch, Chem. –Eur. J., 2018, 24, 7340–7344.
22. P. Cysewski, Phys. Chem. Chem. Phys., 2011, 13, 12998.
23. J. Dominikowska and M. Palusiak, ChemPhysChem, 2018, 19, 590–595.
24. M. Li, X. Wan, X. He, C. Rong and S. Liu, Phys. Chem. Chem. Phys., 2023, 25, 2595–2605.
25. M. C. Caputo and P. Lazzeretti, Int. J. Quantum Chem., 2011, 111, 772–779.
26. T. J. Irons, G. David and A. M. Teale, J. Chem. Theory Comput., 2021, 17, 2166–2185.
27. N. Casati, A. Kleppe, A. P. Jephcoat and P. Macchi, Nat. Commun., 2016, 7, 10901.
28. S. Fanetti, M. Citroni and R. Bini, J. Phys. Chem. C, 2014, 118, 13764–13768.
29. I. I. Naumov and R. J. Hemley, Acc. Chem. Res., 2014, 47, 3551–3559.
30. T. Novoa, J. Contreras-García, P. Fuentealba and C. Cárdenas, J. Chem. Phys., 2019, 150, 204304.
31. X. Dong, A. R. Oganov, H. Cui, X.-F. Zhou and H.-T. Wang, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2117416119.
32. J. Eeckhoudt, T. Bettens, P. Geerlings, R. Cammi, B. Chen, M. Alonso and F. De Proft, Chem. Sci., 2022, 13, 9329–9350.
33. F. Zeller, C.-M. Hsieh, W. Dononelli and T. Neudecker, Wiley Interdiscip. Rev. Comput. Mol. Sci., 2024, 14, e1708.
34. T. Stauch, J. Chem. Phys., 2020, 153, 134503.
35. J. Bentrup, R. Weiß, F. Zeller and T. Neudecker, J. Comput. Chem., 2025, 46, e70024.
36. M. Scheurer, A. Dreuw, E. Epifanovsky, M. Head-Gordon and T. Stauch, J. Chem. Theory Comput., 2020, 17, 583–597.
37. F. Zeller, E. Berquist, E. Epifanovsky and T. Neudecker, J. Chem. Phys., 2022, 157, 184802.
38. A. Pausch, F. Zeller and T. Neudecker, J. Chem. Theory Comput., 2025, 21, 747–761.
39. O. Nielsen and R. M. Martin, Phys. Rev. Lett., 1983, 50, 697.
40. O. Nielsen and R. M. Martin, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 32, 3780.
41. E. Epifanovsky, A. T. Gilbert, X. Feng, J. Lee, Y. Mao, N. Mardirossian, P. Pokhilko, A. F. White, M. P. Coons and A. L. Dempwolff, et al., J. Chem. Phys., 2021, 155, 084801.
42. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
43. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
44. T. H. Dunning Jr, J. Chem. Phys., 1989, 90, 1007–1023.
45. J. Eeckhoudt, M. Alonso, P. Geerlings and F. De Proft, J. Chem. Theory Comput., 2024, 20, 7430–7442.
46. G. Van Rossum, Python 3 Reference Manual, CreateSpace, 2009.
47. J. Pedersen and K. V. Mikkelsen, RSC Adv., 2022, 12, 2830–2842.
48. L. Van Nyvel, M. Alonso and M. Solà, Chem. Sci., 2025, 16, 5613–5622.
49. T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51–57.
50. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
51. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558.
52. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
53. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
54. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
55. J. Kruszewski and T. Krygowski, Tetrahedron Lett., 1972, 13, 3839–3842.
56. T. M. Krygowski, J. Chem. Inf. Comput. Sci., 1993, 33, 70–78.
57. P. v. R. Schleyer, C. Maerker, A. Dransfeld, H. Jiao and N. J. van Eikema Hommes, J. Am. Chem. Soc., 1996, 118, 6317–6318.
58. G. Monaco and R. Zanasi, J. Phys. Chem. A, 2014, 118, 1673–1683.
59. J. Jusélius, D. Sundholm and J. Gauss, J. Chem. Phys., 2004, 121, 3952–3963.
60. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, J. A. Montgomery, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian, J. B. Cross, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, C. Gonzalez and J. A. Pople, Gaussian16 Revision B.03, 2003.
61. R. F. Bader, Atoms in Molecules: A Quantum Theory, Clarendon Press, 1994.
62. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
63. T. A. Keith, TK Gristmill Software: Overland Park, KS, USA, 2019, vol. 23.
64. P. Salvador, E. Ramos-Cordoba, M. Montilla, L. Pujal and M. Gimferrer, J. Chem. Phys., 2024, 160, 172502.
65. E. Matito, ESI-3D: Electron Sharing Indices Program for 3D Molecular Space Partitioning, Institute for Computational Chemistry and Catalysis, University of Girona, 2006.
66. E. Cox, Rev. Mod. Phys., 1958, 30, 159.
67. I. S. Konovalova, S. V. Shishkina, G. Bani-Khaled, E. N. Muzyka and A. N. Boyko, CrystEngComm, 2019, 21, 2908–2919.
68. C. H. Choi and M. Kertesz, J. Chem. Phys., 1998, 108, 6681–6688.
69. M. K. Cyrański and T. M. Krygowski, Tetrahedron, 1999, 55, 6205–6210.
70. S. Shaik, A. Shurki, D. Danovich and P. C. Hiberty, Chem. Rev., 2001, 101, 1501–1540.
71. S. C. Pierrefixe and F. M. Bickelhaupt, Chem. –Eur. J., 2007, 13, 6321–6328.
72. P. v. R. Schleyer and H. Jiao, Pure Appl. Chem., 1996, 68, 209–218.
73. I. F. Perepichka and D. F. Perepichka, Handbook of Thiophene-Based Materials: Applications in Organic Electronics and Photonics, 2 Volume Set, John Wiley & Sons, 2009.
74. T. M. Krygowski and M. K. Cyrański, Chem. Rev., 2001, 101, 1385–1420.
75. J. C. Dobrowolski and S. Ostrowski, J. Chem. Inf. Model., 2023, 63, 7744–7754.
76. E. Matito, J. Poater, M. Duran and M. Solà, J. Mol. Struct. Theochem., 2005, 727, 165–171.
77. E. M. Arpa, S. Stafstrom and B. Durbeej, J. Org. Chem., 2025, 90, 1297–1308.
78. N. C. Craig, P. Groner and D. C. McKean, J. Phys. Chem. A, 2006, 110, 7461–7469.
79. L. Pauling and G. W. Wheland, J. Chem. Phys., 1933, 1, 362–374.
80. M. K. Cyrański, Chem. Rev., 2005, 105, 3773–3811.
81. S. E. Wheeler, K. N. Houk, P. v. R. Schleyer and W. D. Allen, J. Am. Chem. Soc., 2009, 131, 2547–2560.
82. P. v. R. Schleyer and F. Pühlhofer, Org. Lett., 2002, 4, 2873–2876.
83. M. Dewar and H. Schmeising, Tetrahedron, 1959, 5, 166–178.
84. I. Fishtik and R. Datta, J. Phys. Chem. A, 2003, 107, 10471–10476.
85. B. Chen, R. Hoffmann and R. Cammi, Angew. Chem., Int. Ed., 2017, 56, 11126–11142.
86. R. Silverstein and F. Webster, Spectrometry Identification of Organic Compounds, John Wiley and Sons, 6th edn, 1996.
87. R. F. Bader, Chem. Rev., 1991, 91, 893–928.
88. S. Howard and T. Krygowski, Can. J. Chem., 1997, 75, 1174–1181.
89. M. Palusiak and T. M. Krygowski, Chem. –Eur. J., 2007, 13, 7996–8006.
90. I. Casademont-Reig, E. Ramos-Cordoba, M. Torrent-Sucarrat and E. Matito, in Aromaticity, Elsevier, 2021, pp. 235–259.
91. K. Ruedenberg, Rev. Mod. Phys., 1962, 34, 326.
92. R. F. Bader and M. E. Stephens, J. Am. Chem. Soc., 1975, 97, 7391–7399.
93. E. Matito, M. Solà, P. Salvador and M. Duran, Faraday Discuss., 2007, 135, 325–345.
94. R. L. Fulton, J. Phys. Chem., 1993, 97, 7516–7529.
95. J. G. Angyan, M. Loos and I. Mayer, J. Phys. Chem., 1994, 98, 5244–5248.
96. X. Fradera, M. A. Austen and R. F. Bader, J. Phys. Chem. A, 1999, 103, 304–314.
97. J. Grèbol-Tomàs, E. Matito and P. Salvador, Chem. –Eur. J., 2024, 30, e202401282.
98. R. F. Bader, A. Streitwieser, A. Neuhaus, K. E. Laidig and P. Speers, J. Am. Chem. Soc., 1996, 118, 4959–4965.
99. J. Poater, X. Fradera, M. Duran and M. Solà, Chem. –Eur. J., 2003, 9, 400–406.
100. E. Matito, M. Duran and M. Sola, J. Chem. Phys., 2005, 122, 014109.
101. F. Feixas, E. Matito, J. Poater and M. Sola, J. Phys. Chem. A, 2007, 111, 4513–4521.
102. M. Giambiagi, M. S. de Giambiagi and K. C. Mundim, Struct. Chem., 1990, 1, 423–427.
103. M. Giambiagi, M. S. de Giambiagi, C. D. dos Santos Silva and A. P. de Figueiredo, Phys. Chem. Chem. Phys., 2000, 2, 3381–3392.
104. D. M. Esquivel, M. S. De Giambiagi, M. Giambiagi and S. Makler, Chem. Phys. Lett., 1979, 64, 131–138.
105. P. Bultinck, R. Ponec and S. Van Damme, J. Phys. Org. Chem., 2005, 18, 706–718.
106. E. Matito, Phys. Chem. Chem. Phys., 2016, 18, 11839–11846.
107. C. García-Fernandez, E. Sierda, M. Abadia, B. Bugenhagen, M. H. Prosenc, R. Wiesendanger, M. Bazarnik, J. E. Ortega, J. Brede and E. Matito, et al., J. Phys. Chem. C, 2017, 121, 27118–27125.
108. P. Geerlings, F. De Proft and W. Langenaeker, Chem. Rev., 2003, 103, 1793–1874.
109. P. Geerlings, in Conceptual Density Functional Theory, Wiley, 2022, pp. 301–323.
110. N. Sablon, F. D. Proft, P. W. Ayers and P. Geerlings, J. Chem. Theory Comput., 2010, 6, 3671–3680.
111. N. Sablon, F. De Proft and P. Geerlings, J. Phys. Chem. Lett., 2010, 1, 1228–1234.
112. N. Sablon, F. De Proft, M. Solà and P. Geerlings, Phys. Chem. Chem. Phys., 2012, 14, 3960–3967.
113. S. Fias, P. Geerlings, P. Ayers and F. De Proft, Phys. Chem. Chem. Phys., 2013, 15, 2882–2889.
114. S. Fias, Z. Boisdenghien, T. Stuyver, M. Audiffred, G. Merino, P. Geerlings and F. De Proft, J. Phys. Chem. A, 2013, 117, 3556–3560.
115. B. Wang, P. Geerlings, C. Van Alsenoy, F. Heider-Zadeh, P. W. Ayers and F. De Proft, J. Chem. Theory Comput., 2023, 19, 3223–3236.
116. Z. Chen, C. S. Wannere, C. Corminboeuf, R. Puchta and P. v. R. Schleyer, Chem. Rev., 2005, 105, 3842–3888.
117. H. Fallah-Bagher-Shaidaei, C. S. Wannere, C. Corminboeuf, R. Puchta and P. v. R. Schleyer, Org. Lett., 2006, 8, 863–866.
118. A. Stanger, Eur. J. Org Chem., 2020, 2020, 3120–3127.
119. P. Cysewski, Phys. Chem. Chem. Phys., 2011, 13, 12998–13008.
120. P. J. Mayer and H. Ottosson, J. Phys. Org. Chem., 2025, 38, e70000.
121. C. Johnson and F. Bovey, J. Chem. Phys., 1958, 29, 1012–1014.
122. J. A. Pople, J. Chem. Phys., 1956, 24, 1111.
123. J. Waugh and R. Fessenden, J. Am. Chem. Soc., 1957, 79, 846–849.
124. I. Morao and F. P. Cossio, J. Org. Chem., 1999, 64, 1868–1874.
125. F. P. Cossio, I. Morao, H. Jiao and P. v. R. Schleyer, J. Am. Chem. Soc., 1999, 121, 6737–6746.
126. C. Foroutan-Nejad, S. Shahbazian, F. Feixas, P. Rashidi-Ranjbar and M. Sola, J. Comput. Chem., 2011, 32, 2422–2431.
127. S. Pelloni and P. Lazzeretti, J. Phys. Chem. A, 2013, 117, 9083–9092.
128. D. G. Farnum and C. F. Wilcox, J. Am. Chem. Soc., 1967, 89, 5379–5383.
129. J. Elvidge and L. Jackman, J. Chem. Soc., 1961, 859.
130. M. Orozco-Ic, A. Restrepo, A. Muñoz-Castro and G. Merino, J. Chem. Phys., 2019, 151, 014102.
131. I. Roncevic, F. J. Leslie, M. Rossmannek, I. Tavernelli, L. Gross and H. L. Anderson, J. Am. Chem. Soc., 2023, 145, 26962–26972.
132. M. Orozco-Ic and D. Sundholm, Phys. Chem. Chem. Phys., 2023, 25, 12777–12782.
133. S. T. Epstein, J. Chem. Phys., 1973, 58, 1592–1595.
134. J. Juselius, D. Sundholm and J. Gauss, J. Chem. Phys., 2004, 121, 3952–3963.
135. D. Sundholm, H. Fliegl and R. J. Berger, Wiley Interdiscip. Rev. Comput. Mol. Sci., 2016, 6, 639–678.
136. L. Leyva-Parra, R. Pino-Rios, D. Inostroza, M. Solà, M. Alonso and W. Tiznado, Chem. –Eur. J., 2023, 30, e202302415.
137. A. Mohajeri and A. Ashrafi, Chem. Phys. Lett., 2008, 458, 378–383.
138. Z. Badri and C. Foroutan-Nejad, Phys. Chem. Chem. Phys., 2016, 18, 11693–11699.

---

Footnote

† These authors contributed equally to this work.

---