Rationalizing polymorphism with local correlation-based methods: a case study of pnictogen molecular crystals

Abstract

A computational workflow is proposed to quantify and rationalize the relative stability of different structures of molecular crystals using cluster models and quantum chemical methods. The Hartree–Fock plus London Dispersion (HFLD) scheme is used to estimate the lattice energy of molecular crystals in various structural arrangements. The fragment-pairwise Local Energy Decomposition (fp-LED) scheme is then employed to quantify the key intermolecular interactions responsible for the relative stability of different crystal structures. The fp-LED scheme provides also in-depth chemical insights by decomposing each interaction into energy components such as dispersion, electrostatics, and exchange. Notably, this analysis requires only a single interaction energy computation per structure on a suitable cluster model. As a case study, two polymorphs of each of the following are considered: naphthyl-substituted dipnictanes (with As, Sb, and Bi as the pnictogen atom) and tris(thiophen-2-yl)bismuthane. The approach outlined offers high accuracy as well as valuable insights for developing design principles to engineer crystal structures with tailored properties, opening up new avenues in the study of molecular aggregates, potentially impacting diverse fields in materials science and beyond.

---

Introduction

Accurate quantification and analysis of intermolecular interactions in molecular crystals are crucial for comprehending the underlying effects that affect the formation of polymorphs, their stability, and phase transitions. This knowledge is essential for predicting and controlling crystal structures, and for designing materials with tailored properties. In pharmaceuticals, different polymorphs can significantly impact the solubility of a drug as well as its stability and bioavailability.¹ Additionally, in fields like molecular electronics and drug design, the arrangement of molecules within a crystal structure influences the performance and properties of the material.^2,3

To develop structure–property relationships of molecular systems, typically their constituting isolated structural motifs are analyzed. However, conclusions drawn from chemical intuition-based structural analyses can be misleading.^4–6 The interaction strength and nature of a given structural motif can vary across different chemical environments. For instance, what appears to be a classical “π–π interaction” might actually be dominated by dispersion forces.⁷ Similarly, pnictogen–π interactions primarily stem from dispersion, accompanied by a tunable donor–acceptor contribution.^8–10 Hence, computing energetic contributions is essential for properly assessing the nature of the chemical interaction.^11,12 Unfortunately, this task is far from trivial. In particular, for molecular organic crystals, nearly all experimentally-known polymorphs are separated by less than 2.5 kcal mol⁻¹,⁶ necessitating the use of accurate quantum chemical methods. However, these methods become computationally very demanding, if not entirely unfeasible, as the system size grows. Additionally, beyond employing accurate methods, a scheme that accurately quantifies the different physical components of the interaction (electrostatic, exchange, dispersion, etc.) is also essential.

A feasible approach to gain insights into the strength of intermolecular interactions in a molecular crystal is to investigate isolated dimers while neglecting other environmental effects.¹³ For example, density functional theory (DFT) dimer analyses reproduce trends in the interactions between naphthyl-substituted dipnictanes, allowing to explain polymorphism.^14,15 However, it should be kept in mind that for systems like closely packed crystal structures with noticeable charge transfer or polarization effects this approach may fail drastically. For example, for an accurate estimation of the lattice energy of a benzene crystal, it is necessary to compute terms in the many-body expansion (MBE) up to at least the four-body level.¹⁶

Efforts toward describing large molecular systems with increasing accuracy and efficiency constitute a growing active research field. MBE-based methods, such as the method of increments that truncate correlation energy,¹⁷ have significantly extended the limits of computational methods. When combined with embedding/periodic schemes and/or local correlation methods, these approaches enable investigation of not only small models but also a wide range of larger structures, including rare-gas crystals, covalent semiconductors, ionic insulators, diamond, metallic and half-metallic solids, fullerenes, graphite, and polymers.^17–21 Especially periodic local versions of MP2 are a big step forward in accurately investigating large crystal systems.^22–24 When periodic boundary conditions are implemented using the projector-augmented wave method for Coupled Cluster methods, even better accuracies are achieved.²⁵ However, this approach still suffers from the steep scaling of Coupled Cluster methods with system size, restricting their use to solid models such as LiH, noble gases, boron nitride sheets, and the diamond and graphite phases of carbon.^26–28 An alternative approach involves using local variants of CCSD(T), which achieve linear or low-order scaling with system size.^29–31 In particular, the popular Domain-based Local Pair Natural Orbital CCSD(T) method [DLPNO-CCSD(T)]^32–40 allows for the computation of energies and properties of systems as large as entire proteins.⁴¹

To rationalize molecular crystal structures in terms of the underlying intermolecular interactions, in addition to having an accurate and efficient quantum chemical method, a scheme that facilitates the interpretation of the energies is also required. Indeed, Symmetry Adapted Perturbation Theory (SAPT)^42–48 and Energy Decomposition Analysis (EDA)^49–51 methods breakdown interaction energies into physically meaningful components, such as electrostatics, induction, exchange-repulsion, and London dispersion. However, these schemes were primarily developed for the interaction of pairs of fragments. The extension of SAPT to multi-fragment systems marks an active era of research, distinguished by exciting developments.^52–61

On the other hand, the Local Energy Decomposition (LED)^62–65 scheme breaks down the energy of DLPNO-based local correlation methods into several physically meaningful terms for both closed–shell and open–shell molecular systems composed of any given number of fragments. The LED energy terms correlate well with the corresponding SAPT terms in the weak-interaction regime.^66,67 As LED is instrumental not only in weak- but also in strong-interaction regimes, it has found applications in diverse fields.^{8–10,62–76} It is worth emphasizing here that, unlike other fragmentation methods that estimate the energy of the entire system by a series of computations on its isolated fragments and fragment combinations, LED requires only one supermolecular interaction energy computation. In addition, when used in conjunction with its recently introduced fragment pairwise extension fp-LED,⁷⁵ this scheme provides a decomposition of the binding energy (e.g., lattice energies, protein–ligand interactions, etc.) into purely fragment-pairwise contributions.

Hence, LED constitutes a basis for the development of new methods that are accurate and efficient. This led to the semi-empirical “Natural Orbital Tied Constructed Hamiltonian” (NOTCH) method⁷⁷ and the non-empirical “Hartree–Fock plus London Dispersion (HFLD)” method.^78,79 HFLD is a variant of the DLPNO-CCSD(T) method for accurately and efficiently quantifying and analyzing non-covalent interaction (NCI) energies of both closed–hell and open–shell systems. It not only extends the applicability limits of its already efficient parent method significantly but also provides consistently an accuracy between those of CCSD and CCSD(T).^78,79 Therefore, HFLD has found widespread applications in large-scale computations on protein,⁷⁸ DNA,⁸⁰ and crystal⁸¹ structures and on the solute–solvent interactions.⁷⁹

In the present study, the HFLD method is used in conjunction with the fp-LED scheme to investigate the factors affecting the structure and stability of two forms of naphthyl-substituted dipnictane (Pn₂Naph₂, Pn = As, Sb, and Bi)^14,15 crystal structures and of tris(thiophen-2-yl)bismuthane (Bi(C₄H₃S)₃)⁸² crystal structures in the context of polymorphism. With this scheme, a quantification and analysis of all pairwise interactions of the central monomer with the neighboring monomers is provided within one single interaction energy calculation. It is worth noting here that some of these crystal structures were investigated previously using DFT-based dimer calculations which are much more laborious to carry out and yield less accurate results.^15,82

In previous studies,^14,15,82,83 we investigated London dispersion forces in various molecular systems containing group 15 elements as dispersion energy donors, with the aim of assessing their potential in crystal engineering and catalysis. Among the compounds studied, the Bi(C₄H₃S)₃ is an excellent showcase to test our novel theoretical approach for the following reasons: (i) it shows enantiotropic phase transition near the room temperature (ca. 250 K) with low transition energy; (ii) the polymorphism relies on a subtle and reversible change from a dispersion interaction with the π-system and a sulfur atom of one ligand. Last but not least, both polymorphs were successfully characterized by X-ray single-crystal structure analysis.⁸² In the case of Pn₂Naph₂, the crystal structures showed different structural arrangements despite being isovalent.^14,15,83

Thus, these two examples were chosen to illustrate how the proposed scheme can be used to rationalize rather small differences in the crystal structures of related compounds and how it can be used to identify and quantify the interactions that give rise to structural diversity observed in polymorphs.

Structures

Structure of dinaphthyl dipnictane clusters

As₂Naph₂ and Bi₂Naph₂ crystal structures were taken from the literature.^14,15 Despite being isovalent, the observed packing for these two systems differs. We label a molecular arrangement similar to that experimentally found for As₂Naph₂ as “A1”, while we label a molecular arrangement similar to the one observed for Bi₂Naph₂ as “A2”. In this study, we examine the key intermolecular interactions in Pn₂Naph₂ crystal structures with different pnictogen atoms (Pn = As, Sb, Bi) in both A1 and A2 forms to understand the origin of the experimentally observed structural variations.

Initially, the experimentally determined crystal structures for As₂Naph₂ and Bi₂Naph₂ were optimized using periodic boundary conditions at the PBE-D3(BJ) level of theory. Then, 14-molecule (see Fig. 1) and 13-molecule (see Fig. 2) clusters were cut out from the resulting optimized A1 and A2 structures, respectively. These clusters include all monomers directly interacting with the central monomer. To assess the effect of the Pn atom on the stability of the crystal structure, the Pn sites were substituted with Sb and Bi in the A1 cluster (originally containing As), while they were substituted with As and Sb in the A2 cluster (originally containing Bi), without any further geometry optimization.

Fig. 1 (a) Pn₂Naph₂ (Pn = As, Sb, and Bi) monomer. (b) Schematic representation of its A1 form crystal together with the labeling of the constituting monomers, which is formed with Pn = As. (c) The corresponding 14-molecule cluster model of the A1 form cut from the PBE-D3(BJ)-optimized crystal structure (vide supra). Color code: H, white; Pn, green; C, pink for the central monomer and gray for the environmental monomers.

Fig. 2 (a) Pn₂Naph₂ (Pn = As, Sb, and Bi) monomer. (b) Schematic representation of its A2 form crystal together with the labeling of the constituting monomers, which is formed with Pn = Bi. (c) The corresponding 13-molecule model of the A2 form cut from the PBE-D3BJ-optimized crystal structure (vide supra). Color code: H, white; Pn, green; C, pink for the central monomer and gray for the environmental monomers.

Structure of tris(thiophen-2-yl)bismuthane clusters

Bi(C₄H₃S)₃ crystallizes in the space group of R3̄ at 269 K (labeled as B1 in this study) and in P1̄ at 245 K (labeled as B2 in this study). These crystal structures show an enantiotropic phase transition at 250 K with a transition energy of 1.4 kJ mol⁻¹.⁸² Although both clusters look very similar, B2 packs more closely than B1.⁸² 17-Molecule clusters were cut out from the supramolecular arrangements as determined by single crystal X-ray structure analyses⁸² of both B1 and B2 polymorphs without any geometry optimization. The resulting clusters that include all monomers directly interacting with the central monomer are shown in Fig. 3 together with the labeling of their monomers.

Fig. 3 (a) Bi(C₄H₃S)₃ monomer. (b) Schematic representation of its crystal structure together with the labeling of its constituting monomers. (c) The corresponding 17-molecule cluster of the B1 form. (d) The corresponding 17-molecule cluster of the B2 form. Color code: H, white; Pn, green; S, yellow; C, pink for the central monomer and gray for the environmental monomers.

Computational scheme and details

Combined HFLD/fp-LED scheme for interaction energies

Labeling the energy of the entire cluster as E, the energy of the central monomer as E₀, and the energy of the rest of the cluster (the “environment”) as E_env, the interaction energy ΔE of the central monomer with the environment becomes within the supermolecular approach:

ΔE = E − (E₀ + E_env) (1)

It is important to note that, in the framework of MBE for computing the lattice energy, ΔE is expressed as the sum of the contributions from all isolated dimers, trimers, tetramers, and so on. When this expansion is truncated at the dimer level for highly symmetric molecular crystals, half of ΔE corresponds to the electronic component of the lattice energy.⁸⁴ This estimate of lattice energy can be systematically refined by incorporating higher-body terms using well-established scaling protocols.⁸⁴ In our approach, rather than summing the contributions from individual fragment combinations, we decompose ΔE computed using the supermolecular approach into pairwise contributions, which inherently include many-body effects. For the crystals considered in this study, by symmetry, ΔE represents twice the electronic component of the lattice energy.

In eqn (1), the entire system and the environment are both composed of multiple monomers, and thus require performing LED analyses on top of their standard electronic structure calculations to obtain pairwise contributions. Using the notation introduced in ref. 75 for the standard LED decomposition of interaction energies between subsystems of many fragments, we obtain:

The summations run over the individual fragments (the monomers) labeled as X, Y = 0, 1, 2,…, (N − 1), where N is the number of fragments in the cluster. The electronic preparation term ΔE^el-prep_X corresponds to the change in the energy of monomer X upon the interaction between the central monomer and the environment. It is shown by the diagonal elements (i.e., X = Y) enclosed by violet boxes in Fig. 4(a).

Fig. 4 A1 form of As₂Naph₂: HFLD/LED interaction energy map of the central monomer 0 with those at its environment (1–13), using (a) standard LED and (b) fp-LED. For simplicity, interaction terms with absolute values smaller than 0.2 kcal mol⁻¹ are not annotated on the heat maps. See ESI† for fully annotated heat maps.

The interaction between fragment X and fragment Y (where X ≠ Y) is denoted as ε^(X,Y), or equivalently as ε^(Y,X). To avoid double counting of ε^(X,Y) and ε^(Y,X), the sum over X and Y was constrained in the second term of eqn (2) with the X > Y condition. In the heat maps, we place X labels on the horizontal axis, while placing Y labels on the vertical axis. This choice corresponds to upper-diagonal heat maps as given in Fig. 4.

The pairwise ε^(X,Y) terms (where X ≠ Y) are of two kinds: genuine interaction terms between the central monomer and each fragment in the environment, ε^(X,0) (elements enclosed by black boxes in the first row of Fig. 4(a)), and “inductive” interaction terms, ε^(X,Y≠0). The inductive terms (elements enclosed by green boxes in Fig. 4(a)) represent the changes in fragment-pairwise interaction energies within the environment that occur upon the formation of the entire cluster. These terms arise from the cooperativity effects of various noncovalent interactions and thus vanish in the absence of many-body effects.

In eqn (2), ΔE^el-prep_X incorporates the dominant repulsive contributions for the interaction of fragment X with all the remaining fragments in the system. The fp-LED⁷⁵ scheme allows to disentangle the contributions from pairs of fragments within ΔE^el-prep_X terms as

ΔE^el-prep_XY = ω^(X,Y)_XΔE^el-prep_X + ω^(X,Y)_YΔE^el-prep_Y (3)

wherewith K = 0, 1, 2,…, (N − 1). The denominators of ω^(X,Y)_X and ω^(X,Y)_Y in eqn (4) are simply the sum of all pairwise energy terms involving fragment X and fragment Y, respectively. These denominators ensure that the sum of all ΔE^el-prep_X terms is equal to the sum of all pairwise ΔE^el-prep_XY terms. Plugging eqn (3) into eqn (2), one obtains fully pairwise decomposed interaction energy:

The standard LED map based on eqn (2) contains positive and negative values of large magnitude (Fig. 4(a)). In contrast, the fp-LED map based on eqn (5) (Fig. 4(b)) directly provides the strength of all interactions of the fragment pairs within the structure, with interaction energy values comparable to those of the isolated fragment pairs. However, unlike the isolated dimer computations, these values also include the effect of the chemical environment, thus provide an exact decomposition of the interaction energy into dimer contributions.

Eqn (3) was formulated based on the perfect correlation between ε^(X,Y) and the sum of ΔE^el-prep_X and ΔE^el-prep_Y for noncovalently interacting isolated X⋯Y dimers.⁷⁵ In other words, the strength of pairwise interaction is proportional to the increase in the energy of the fragments upon their interaction. This proportionality suggests that, by design, the ε^(X,Y) and ΔE_int,XY terms in standard and fp-LED typically follow the same trends. This is also reflected in the LED data for all clusters considered in this study (see the ESI†). Hence, despite the fact that interpreting standard LED results is less intuitive due to very large ε^(X,Y) and ΔE^el-prep_X values, which often have opposite signs, standard LED is still useful in trend studies.

In the HFLD scheme,^78–80 the ΔE_int,XY terms contain electronic preparation, electrostatic, and exchange components computed at the HF level as well as the dispersion component from Coupled Cluster. In this study, for the sake of simplicity, ΔE_int,XY terms are simply decomposed into non-dispersive and dispersive components, which we denote as ΔE_HF,XY and E_disp,XY. Decomposed nondispersive pairwise terms (electronic preparation, electrostatic, and exchange) are all provided in the ESI.† In the HFLD scheme, inductive E_disp,XY dispersion terms (where Y ≠ 0 in the present study) are neglected, as they are practically zero for noncovalent interactions.⁸⁰ This significantly enhances the efficiency of the method without sacrificing accuracy. Therefore, only genuine E_disp,XY dispersion terms (where Y = 0 in the present study) are computed within the HFLD scheme.

It is also worth mentioning that the “inductive” ε^(X,Y≠0) terms, depicted in green in Fig. 4(b), are very small, summing to less than 1 kcal mol⁻¹ in magnitude. Therefore, in the following, we will only display the first row of the LED maps, plus a cell labeled “rest” at the end of the maps that corresponds to the sum of all of these cooperative ε^(X,Y≠0) terms. Full LED maps, as shown in Fig. 4, are provided in the ESI.†

Finally, it is important to emphasize that the pairwise terms just discussed are not derived from isolated dimer calculations but are obtained through a single supermolecular interaction energy calculation, and thus account for many-body effects. The supermolecular approach involves electronic structure calculations only on the central molecule and two molecular clusters: the environment surrounding the central molecule and the entire system. Note that the HFLD scheme^78,79 for noncovalent interactions already efficiently handles molecular clusters of ∼1000 atoms with triple-ζ quality basis sets (∼14 000 basis functions).^74,75,80

Computational details

The experimentally determined crystal structures for As₂Naph₂¹⁴ and Bi₂Naph₂¹⁵ were initially optimized in the solid-state with periodic boundary conditions using the PWSCF (v.6.7) module of the QuantumESPRESSO package.^85,86 The PBE⁸⁷ functional was used by incorporating the D3⁸⁸ dispersion correction with Becke–Johnson damping (BJ).⁸⁹ Ultrasoft Vanderbilt-type pseudo potentials were utilized.⁹⁰ A kinetic-energy cutoff of 40 Ry and a 7 × 4 × 2 Monkhorst–Pack k-point grid⁹¹ was employed.

All the multi-fragment HFLD/LED computations were performed with the ORCA program package (version 5.0)^92–95 at the cluster geometries described in the Structures section. The RIJCOSX approach^96,97 was used for the reference HF part. The def2-TZVP(-f) basis set was used together with its matching auxiliary basis sets in the HF and correlation parts.⁹⁸ NormalPNO* settings⁷⁸ of DLPNO Coupled Cluster methods and default frozen core settings⁹⁹ were applied for the correlation energy calculations. The Foster-Boys scheme¹⁰⁰ was used for localizing both occupied orbitals and PNOs.

Since in the HFLD scheme the electronic preparation is computed at the HF level, the ω^(X,Y)_X and ω^(X,Y)_Y terms in eqn (3) were also computed at the HF level (neglecting the dispersion contribution), as suggested in the theory paper of fp-LED.⁷⁵ LED Analysis Wizard (LEDAW) program package¹⁰¹ was used in computing interaction energy matrices from ORCA outputs and plotting the corresponding heat maps.

Dispersion Interaction Density (DID)^64,102 plots were computed to assess the spatial origin of the dispersion interaction between the central monomer and environmental monomers by using the above-described HFLD/LED settings. For the ease of visualization, DID plots are shown herein for the monomer pairs that most significantly contribute to the stability of the crystal structures separately rather than on the entire cluster.

Results and discussion

In this section, on two forms of Pn₂Naph₂ (A1 and A2) and Bi(C₄H₃S)₃ (B1 and B2) crystals we assess (i) the energetic contribution of each dimeric motif on the entire cluster formation, (ii) the effect of the type of Pn atom on the individual noncovalent interactions and on the overall stability of a given cluster; (iii) the individual interactions responsible for the difference in the binding energies, and thus, lattice energies of the two forms.

Dinaphthyl dipnictane clusters

Pn₂Naph₂ molecules form crystal structures with quite diverse intermolecular interaction patterns.^14,15,83 As₂Naph₂¹⁴ (see the A1 form in Fig. 1) exhibits an intermolecular network of π⋯π, CH⋯π, As⋯As interactions,¹⁵ whereas Bi₂Naph₂ (see the A2 form in Fig. 2) features two Bi⋯π contacts per Bi atom rather than π⋯π contacts in the solid state.¹⁵ Through a very large series of DFT-D dimer computations,¹⁵ the A2 model with Pn = Bi was shown to benefit from the enhanced Pn⋯π intermolecular dispersion interaction. In order to demonstrate the advantage of combined HFLD/fp-LED scheme that allows quantifying and rationalizing all intermolecular interactions in crystal structures with one supermolecular interaction energy calculation, in this study we revisit these A1 and A2 models together with a series of hypothetical structures obtained by just replacing the Pn atom by their lighter or heavier homologs (see the Structures section).

The HFLD/fp-LED pairwise interaction energy map and its nondispersive HF and dispersive interaction energy sub-maps are shown for the A1 form (Pn = As) and its Sb- and Bi-substituted variants in Fig. 5 and for the A2 form (Pn = Bi) and its As- and Sb-substituted variants in Fig. 6. The total HFLD interaction energy ΔE as well as its total nondispersive HF interaction energy ΔE_HF and the total dispersion energy E_disp components are provided in panel b of these figures while their pairwise components are provided in panel c–e.

Fig. 5 (a) Schematic representation of the A1 polymorph of Pn₂Naph₂ with the labeling of their monomers; (b) HFLD interaction energy and its nondispersive HF and dispersion components; HFLD/fp-LED interaction energy maps of the central (c) As₂Naph₂, (d) Sb₂Naph₂, and (e) Bi₂Naph₂ monomer with those in their environment together with nondispersive (HF) and dispersive interaction energy sub-maps. All energy values are given in kcal mol⁻¹. Red denotes attractive and blue denotes repulsive interaction.

Fig. 6 (a) Schematic representation of the A2 polymorph of Pn₂Naph₂ with the labeling of their monomers; (b) HFLD interaction energy and its nondispersive HF and dispersion components; HFLD/fp-LED interaction energy maps of the central (c) As₂Naph₂, (d) Sb₂Naph₂, and (e) Bi₂Naph₂ monomer of these structures with those at their environment together with nondispersive (HF) and dispersive interaction energy sub-maps. All energy values are given in kcal mol⁻¹. Red denotes attractive and blue denotes repulsive interaction.

For the A1 form, with increasing size of the Pn atom, overall repulsive non-dispersive and overall attractive dispersion energies both increase roughly at the same rate. This trend is due to the π⋯π interaction between one of the naphthyl moieties of the central monomer 0 and monomer 8 (and, to a smaller extent, monomer 6). The other pairwise interaction energies remain roughly the same with the change of the Pn atom. Consequently, the sum of all interaction energies does not exhibit a strong dependence on the type of the Pn atom for the A1 form.

The primary stabilizing pairwise interaction in the A1 form emerges from the π⋯π (appearing once in the pair) and CH⋯π (appearing once in the pair) interactions between one of the naphthyl moieties of the central monomer 0 and the monomer 4, with a secondary contribution from the same type of interaction between central monomer 0 and monomer 6 (see Fig. 7(a) for the structural motif of these monomer pairs), consistent with previous DFT-D results.¹⁵ These interactions are mainly of dispersion nature and the associated DID plot is shown in Fig. 7(a). The combined contribution of these two pairs constitutes ca. 50% of the total interaction energy. Although individual contributions from other pairs are much smaller, their collective effect is significant, accounting for the remaining half of the crystal structure's stability.

Fig. 7 The structural motif that has the largest contribution to the stability of the polymorphs of Pn₂Naph₂ and its DID plot at the HFLD/LED level (a) on the A1 form that features π⋯π and CH⋯π interactions (b) on the A2 form that features Bi⋯π interactions. These structural motifs appear twice in the cluster models of the crystal structures. The DID isosurface value was set to 0.004 kcal mol⁻¹ bohr⁻³.

Transition from A1 to A2 involves a decrease in both the nondispersive HF and dispersive interaction components, with exception of the dispersion component in the case that Pn = Bi. Therefore, the enhanced stability of the A2 form with Pn = Bi compared with A1 stems from a dual effect: a reduction in the nondispersive component (ca. 50%) and an increase in the dispersive interaction (ca. 50%).

Unlike A1, the overall stability of A2 significantly increases with the increase in the size of the Pn atom, governed primarily by the increase in the dispersion interaction. This is actually consistent with the larger dispersion energy donation capacity of the larger Pn atoms.

The interaction energy of A2 is only marginally larger than that of A1 when considering Pn = As, yet it demonstrates significantly greater stability with Pn = Sb and Bi. This implies the potential for synthesizing the A2 form also with Pn = As and Sb under diverse experimental conditions. Notably, in the case of Pn = Bi, consistent with the experimental isolation of the Bi structure in the A2 form, the stability of the A2 crystal structure surpasses that of A1 by approximately 20 kcal mol⁻¹.

The pairs that align precisely on top of each other, i.e., the 0–4 and 0–6 monomer pairs, provide the largest pairwise contribution to the stability of the A2 form (see Fig. 7(b) for their structural motif), consistent with the results of the previous DFT-D study.¹⁵ These pairs exhibit Pn⋯π interactions between the Pn atom of one monomer and both naphthyl moieties of the other monomer. This interaction is of mainly dispersive character and increases with the size of the Pn atom (see Fig. 7(b) for the associated DID plot). These two fragment pairs collectively contribute 50%, 52%, and 57% to the total interaction energy with Pn = As, Sb, and Bi, respectively. While the contributions from other individual fragment pairs are much smaller, their cumulative impact becomes notably significant.

As a closing remark of this section, it is worth emphasizing that the formation of both polymorphs, A1 and A2, is primarily driven by dispersive forces irrespective of the type of Pn atom as their nondispersive HF contributions are repulsive. Hence, using the approach described above, it is possible to quantify the relative stability of a given crystal structure, to identify individual structure determining interactions and explain their contribution using the concepts of intermolecular interactions.

Tris(thiophen-2-yl)bismuthane clusters

Two interconvertible polymorphs of Bi(C₄H₃S)₃ were previously crystallized at distinct temperatures.⁸² We label here the high temperature polymorph as B1 while labeling the low temperature polymorph as B2 (see Fig. 3 for their structures and monomer labeling).

The polymorphs show very small but obviously decisive structural changes in the arrangement of the molecules in the unit cell. Based on our previous DFT-D and MP2 studies on small model systems,⁸² the Bi⋯π and Bi⋯S interactions are primarily dispersive in nature. Therefore, the formation of the Bi(C₄H₃S)₃ structure should also be governed primarily by dispersion forces. Furthermore, gas-phase geometry optimizations on the Bi(C₄H₃S)₃⋯(C₄H₃S)_n system (n = 1 and 3) yielded only one minimum.⁸² This suggests that the polymorphism of the Bi(C₄H₃S)₃ system stems from the interplay between different intermolecular interactions, which induce several structural rearrangements and give rise to polymorphism.

In order to shed light into the origin of the polymorphism, we performed HFLD/fp-LED computations on the large Bi(C₄H₃S)₃ clusters of B1 and B2 (see Fig. 3). The resulting HFLD/fp-LED pairwise interaction energy map and its nondispersive HF and dispersive interaction energy sub-maps are shown in Fig. 8 together with the total HFLD interaction energy ΔE, its total nondispersive HF interaction energy ΔE_HF, and total dispersion energy E_disp components.

Fig. 8 (a) Schematic representation of the B1 and B2 polymorphs of Bi(C₄H₃S)₃ with the labeling of their monomers; (b) HFLD interaction energy of the structures and its nondispersive HF and dispersion components; HFLD/fp-LED interaction energy maps of the central monomer of the (c) B1 and (d) B2 forms of Bi(C₄H₃S)₃ with those at its environment together with nondispersive (HF) and dispersive interaction energy sub-maps. All energy values are given in kcal mol⁻¹. Red denotes attractive and blue denotes repulsive interaction.

The B2 cluster is 4.6 kcal mol⁻¹ more stable than the B1 cluster. The overall nondispersive HF interaction components of both B1 and B2 are repulsive and 3–4 times smaller in absolute value than their overall attractive dispersion components. Therefore, the formation of these structures is largely driven by dispersion interaction between Bi(C₄H₃S)₃ molecules, consistent with the results on small model systems.⁸²

The largest contribution to the stability of these polymorphs arises from the interaction of the central monomer 0 with monomers 13, 9, 7, 8, and 1, listed in descending interaction strength on B2. These pairs collectively account for 69% and 72% of the total interaction energy of B1 and B2. Although the cumulative impact of the remaining 11 pairs to the stability of the clusters is still significant, their individual contributions do not vary much between B1 and B2. Therefore, in the following we focus on the most strongly interacting five monomer pairs listed above.

To assess the relation between the structure and the stability of the monomer pairs, in Fig. 9(a) we show the superimposed structures (with respect to the central monomer) of the most stable five monomer pairs in the B1 and B2 polymorphs. As the stabilizing component of the interaction is largely intermolecular dispersion, we also provide the DID plots of these pairs on the B1 and B1 clusters in Fig. 9(b) and (c), respectively.

Fig. 9 (a) The most contributing structural motifs to the stability of the B1 and B2 polymorphs of Bi(C₄H₃S)₃ crystal structures superimposed by minimizing RMS deviation between the Cartesian coordinates of their central monomers indexed 0. For the ease of distinction, C atoms in B1 and B2 are shown pink and gray, respectively. (b) For B1 and (c) for B2 monomer pairs, DID plots computed at the HFLD/LED level. The DID isosurface value was set to 0.003 kcal mol⁻¹ bohr⁻³.

Transitioning from B1 to B2, the distance between the monomers of the strongest pair (0–13) is significantly shortened, including the intermolecular C₄H₃S⋯C₄H₃S ring and the Bi⋯Bi pnictogen separations. This introduces a slight increase in the repulsive component of the interaction while notably enhancing dispersive CH⋯π interactions between the rings (see Fig. 8 for the energy contributions and Fig. 9 for the associated DID plots). These structural modifications account for more than half (3.0 kcal mol⁻¹) of the larger stability of B2 compared to B1.

Pairs 0–7, 0–8, and 0–9, involving intermolecular π⋯π, Bi⋯π, Bi⋯S, π⋯S, S⋯S interactions (see Fig. 9), experience noticeably large repulsion, which accounts for 79% and 72% (see Fig. 8) of the overall repulsive contribution to the interaction energy of B1 and B2, respectively. Despite this, the magnitude of the dispersion interaction in these pairs is akin to that computed for the strongest pair (0–13) of B2. Consequently, due to the counteracting effects of increased and decreased attractive and repulsive components, these pairs demonstrate individual interaction strengths similar to those in the 0–13 pair of B1 discussed above. While transitioning from B1 to B2, as the structural variations in these pairs are relatively minor, they do not noticeably contribute to the relative stability of B1 and B2, despite some variations in their individual interaction strengths and their components.

The interaction strength of the 0–1 pair is significantly smaller than the other four pairs discussed above. During the transition from B1 to B2, the intermolecular distance between its monomers is noticeably shortened. This enhances intermolecular dispersive CH⋯π and S⋯S interactions (see Fig. 9) between the monomers by 2.3 kcal mol⁻¹ (see Fig. 8), accounting for the other half of the larger stability of B2.

To sum up, a combined HFLD/fp-LED analysis allows to quantify and analyze the energetic contributions to the stability of different polymorphs in Bi(C₄H₃S)₃ crystal structures. The analysis demonstrates that the larger stability of the polymorph B2 over B1 arises from the increase of the (predominantly dispersion) interaction of central monomer 0 with monomers 13 and 1, facilitated by their relatively shorter intermolecular contacts.

Conclusions

In this study, HFLD/fp-LED computations were conducted on large Pn₂Naph₂ and Bi(C₄H₃S)₃ cluster models of crystal structures to estimate the binding energy of the central monomer in various possible polymorphs, which are related to lattice energies. These computations provide insights into the relative stability of the polymorphs. The proposed protocol provides an exact decomposition of approximated binding energies into contributions from pairs of monomers within the solid. This computational approach enabled us to rationalize the relative stability and polymorphism of the structures based on fundamental energetic considerations, rather than relying solely on the analysis of their geometry parameters, which are not necessarily correlated with stability.

For the Pn₂Naph₂ polymorphs (A1 and A2), our present HFLD/fp-LED binding energy computations revealed that their most stable dimer structure accounts for half of the total stability of each polymorph, consistent with our earlier DFT-D results on isolated dimers.¹⁵ While the individual contributions of other pairs are much smaller, their cumulative effect accounts for the remaining stability. The enhanced stability of the A2 polymorph with Pn = Bi compared to the A1 polymorph is consistent with the experimental isolation of only A2 with Bi.¹⁵ This stabilization arises equally from a reduction in the nondispersive component and from an increase in the dispersive component. In the case of Pn = As, although only A1 has been experimentally isolated,¹⁴A1 and A2 are computed to be nearly isoenergetic. This suggests the potential for preparing A2 under different experimental conditions.

Having two polymorphs of Bi(C₄H₃S)₃ (B1 and B2) at hand, we were able to gather insights into an enantiotropic phase transition having a low transition barrier. Our HFLD/fp-LED approach allowed us, for the first time, to decipher the subtle interplay of different intermolecular interactions on the stability of each polymorph and their relative strengths. In particular, we demonstrated that subtle changes in the relative orientation of two dimers in these polymorphs have a strong influence on their dispersion interactions, and are therefore almost entirely responsible for their relative stability.

In conclusion, the combined HFLD/fp-LED scheme is a powerful tool for analyzing structural motifs in molecular crystals by decomposing interaction energies, and thus lattice energies, into pairwise intermolecular interaction components. This approach offers valuable insights into the key factors that influence the formation and stability of molecular crystals and their polymorphs. Understanding how molecules aggregate into diverse structures is a significant advancement, especially considering that experiments often yield unexpected structures with slight variations in crystallization conditions. By rationalizing the effects of subtle changes in geometry and electronic structure on crystal stability through the HFLD/fp-LED scheme, material design strategies can be significantly enhanced. Given the pivotal role of crystal structure in various fields, from molecular electronics to drug design, this knowledge holds great potential for driving innovation across various industries.

Finally, it is worth noting that ongoing efforts in our theory lab to extend the LED scheme to covalent interactions are expected to significantly broaden the applicability of this approach to a wider range of molecular systems and interactions.

Author contributions

E. S. and A. A. carried out all the calculations. A. A. analyzed the results and wrote the original draft. A. A. A. and G. B. devised and coordinated the project, and prepared the final version of the manuscript with input from all authors. S. S. and M. M. were involved in the project design and writing and proof reading of the manuscript.

Data availability

The XYZ coordinates and energetic data of all clusters, including heat maps, are available in the ESI.†

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors gratefully acknowledge the Priority Program “Control of Dispersion Interactions in Molecular Chemistry” (SPP 1807) of the Deutsche Forschungsgemeinschaft for financial support. A. A. gratefully acknowledges the Research Unit Program “Bioinspired Oxidation Catalysis with Iron Complexes” (FOR 5215) of the Deutsche Forschungsgemeinschaft for financial support. G. B. acknowledge financial support under the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.1, Call for tender no. 104 published on 2.2.2022 by the Italian Ministry of University and Research (MUR), funded by the European Union – NextGenerationEU – Project Title 2022WKTH9E – Modelling and design of organic conjugated redox materials for energy-saving applications: a bottom-up strategy. – CUP J53D23008810006 – Grant Assignment Decree No. 1064 adopted on 18.07.2022 by the Italian Ministry of University and Research (MUR). Open Access funding provided by the Max Planck Society.

References

1. J. Bernstein, Polymorphism in Molecular Crystals, Oxford University Press, Oxford, 2nd edn, 2020.
2. H. Chung and Y. Diao, J. Mater. Chem. C, 2016, 4, 3915–3933.
3. M. Li, A. H. Balawi, P. J. Leenaers, L. Ning, G. H. L. Heintges, T. Marszalek, W. Pisula, M. M. Wienk, S. C. J. Meskers, Y. Yi, F. Laquai and R. A. J. Janssen, Nat. Commun., 2019, 10, 1–11.
4. B. K. Saha, N. K. Nath and R. Thakuria, Chem. Rec., 2023, 23, e202200173.
5. S. Varughese, J. Mater. Chem. C, 2014, 2, 3499–3516.
6. G. J. O. Beran, Chem. Sci., 2023, 14, 13290–13312.
7. S. Grimme, Angew. Chem., Int. Ed., 2008, 47, 3430–3434.
8. E. Schiavo, K. Bhattacharyya, M. Mehring and A. A. Auer, Chem. – Eur. J., 2021, 27, 14520–14526.
9. M. Krasowska, W. B. Schneider, M. Mehring and A. A. Auer, Chem. – Eur. J., 2018, 24, 10238–10245.
10. M. Krasowska, A. M. Fritzsche, M. Mehring and A. A. Auer, ChemPhysChem, 2019, 20, 2539–2552.
11. S. Rindfleisch, M. Krull, J. Uranga, T. Schmidt, F. Rabe von Pappenheim, L. L. Kirck, A. Balouri, T. Schneider, A. Chari, R. Kluger, G. Bourenkov, U. Diederichsen, R. A. Mata and K. Tittmann, Nat. Catal., 2022, 5, 332–341.
12. X. Aniban, B. Hartwig, A. Wuttke and R. A. Mata, Phys. Chem. Chem. Phys., 2021, 23, 12093–12104.
13. A. M. Preda, M. Krasowska, L. Wrobel, P. Kitschke, P. C. Andrews, J. G. MacLellan, L. Mertens, M. Korb, T. Rüffer, H. Lang, A. A. Auer and M. Mehring, Beilstein J. Org. Chem., 2018, 14, 2125–2145.
14. K. Dzialkowski, A. Gehlhaar, C. Wölper, A. A. Auer and S. Schulz, Organometallics, 2019, 38, 2927–2942.
15. A. Gehlhaar, E. Schiavo, C. Wölper, Y. Schulte, A. A. Auer and S. Schulz, Dalton Trans., 2022, 51, 5016–5023.
16. J. Yang, W. Hu, D. Usvyat, D. Matthews, M. Schütz and G. K. L. Chan, Science, 2014, 345, 640–643.
17. H. Stoll, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 6700.
18. E. Voloshina, N. Gaston and B. Paulus, J. Chem. Phys., 2007, 126, 134115.
19. B. Paulus, Phys. Rep., 2006, 428, 1–52.
20. U. Birkenheuer, P. Fulde and H. Stoll, Theor. Chem. Acc., 2006, 116, 398–403.
21. H. Stoll, B. Paulus and P. Fulde, J. Chem. Phys., 2005, 123, 144108.
22. L. Maschio, B. Civalleri, P. Ugliengo and A. Gavezzotti, J. Phys. Chem. A, 2011, 115, 11179–11186.
23. L. Maschio, D. Usvyat and B. Civalleri, CrystEngComm, 2010, 12, 2429–2435.
24. C. Pisani, M. Schütz, S. Casassa, D. Usvyat, L. Maschio, M. Lorenz and A. Erba, Phys. Chem. Chem. Phys., 2012, 14, 7615–7628.
25. I. Y. Zhang and A. Grüneis, Front. Mater., 2019, 6, 432749.
26. S. J. Nolan, M. J. Gillan, D. Alfè, N. L. Allan and F. R. Manby, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 165109.
27. T. Gruber, K. Liao, T. Tsatsoulis, F. Hummel and A. Grüneis, Phys. Rev. X, 2018, 8, 021043.
28. A. Grüneis, G. H. Booth, M. Marsman, J. Spencer, A. Alavi and G. Kresse, J. Chem. Theory Comput., 2011, 7, 2780–2785.
29. Linear-Scaling Techniques in Computational Chemistry and Physics: Methods and Applications, ed. R. Zaleśny, M. G. Papadopoulos, P. G. Mezey and J. Leszczynski, Springer, 2011.
30. Q. Ma, M. Schwilk, C. Köppl and H.-J. Werner, J. Chem. Theory Comput., 2017, 13, 4871–4896.
31. M. Nagy and P. R. Kállay, J. Chem. Theory Comput., 2019, 15, 5275–5298.
32. F. Neese, A. Hansen and D. G. Liakos, J. Chem. Phys., 2009, 131, 064103.
33. F. Neese, A. Hansen, F. Wennmohs and S. Grimme, Acc. Chem. Res., 2009, 42, 641–648.
34. F. Neese, F. Wennmohs and A. Hansen, J. Chem. Phys., 2009, 130, 114108.
35. A. Hansen, D. G. Liakos and F. Neese, J. Chem. Phys., 2011, 135, 214102.
36. L. M. J. Huntington, A. Hansen, F. Neese and M. Nooijen, J. Chem. Phys., 2012, 136, 064101.
37. C. Riplinger and F. Neese, J. Chem. Phys., 2013, 138, 034106.
38. D. G. Liakos, A. Hansen and F. Neese, J. Chem. Theory Comput., 2011, 7, 76–87.
39. C. Riplinger, P. Pinski, U. Becker, E. F. Valeev and F. Neese, J. Chem. Phys., 2016, 144, 024109.
40. D. G. Liakos, M. Sparta, M. K. Kesharwani, J. M. L. Martin and F. Neese, J. Chem. Theory Comput., 2015, 11, 1525–1539.
41. C. Riplinger, B. Sandhoefer, A. Hansen and F. Neese, J. Chem. Phys., 2013, 139, 134101.
42. B. Jeziorski, R. Moszynski and K. Szalewicz, Chem. Rev., 1994, 94, 1887–1930.
43. R. Podeszwa, R. Bukowski and K. Szalewicz, J. Chem. Theory Comput., 2006, 2, 400–412.
44. A. Heßelmann, G. Jansen and M. Schütz, J. Chem. Phys., 2005, 122, 014103.
45. R. Bukowski, R. Podeszwa and K. Szalewicz, Chem. Phys. Lett., 2005, 414, 111–116.
46. A. Heßelmann, J. Chem. Theory Comput., 2018, 14, 1943–1959.
47. B. Jeziorski, M. van Hemert and B. Jeziorski, Mol. Phys., 1976, 31, 713–729.
48. T. M. Parker, L. A. Burns, R. M. Parrish, A. G. Ryno and C. D. Sherrill, J. Chem. Phys., 2014, 140, 094106.
49. E. Pastorczak and C. Corminboeuf, J. Chem. Phys., 2017, 146, 120901.
50. M. von Hopffgarten and G. Frenking, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 43–62.
51. M. J. S. Phipps, T. Fox, C. S. Tautermann and C.-K. Skylaris, Chem. Soc. Rev., 2015, 44, 3177–3211.
52. R. M. Parrish, K. C. Thompson and T. J. Martínez, J. Chem. Theory Comput., 2018, 14, 1737–1753.
53. R. M. Parrish, T. M. Parker and C. David Sherrill, J. Chem. Theory Comput., 2014, 10, 4417–4431.
54. R. M. Parrish and C. D. Sherrill, J. Chem. Phys., 2014, 141, 44115.
55. A. Heßelmann, Beilstein J. Org. Chem., 2018, 14, 979–991.
56. J. Garcia, R. Podeszwa and K. Szalewicz, J. Chem. Phys., 2020, 152, 184109.
57. D. Luu, C. Corminboeuf and K. Patkowski, J. Chem. Theory Comput., 2024, 20, 7884–7903.
58. D. Luu and K. Patkowski, J. Phys. Chem. A, 2023, 127, 356–377.
59. O. R. Meitei and A. Heßelmann, J. Comput. Chem., 2017, 38, 2500–2508.
60. E. Pastorczak, A. Prlj, J. F. Gonthier and C. Corminboeuf, J. Chem. Phys., 2015, 143, 224107.
61. J. F. Gonthier and C. Corminboeuf, J. Chem. Phys., 2014, 140, 154107.
62. W. B. Schneider, G. Bistoni, M. Sparta, M. Saitow, C. Riplinger, A. A. Auer and F. Neese, J. Chem. Theory Comput., 2016, 12, 4778–4792.
63. A. Altun, M. Saitow, F. Neese and G. Bistoni, J. Chem. Theory Comput., 2019, 15, 1616–1632.
64. A. Altun, F. Neese and G. Bistoni, J. Chem. Theory Comput., 2019, 15, 215–228.
65. G. Bistoni, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2019, e1442.
66. A. Altun, F. Neese and G. Bistoni, Beilstein J. Org. Chem., 2018, 14, 919–929.
67. A. Altun, R. Izsák and G. Bistoni, Int. J. Quantum Chem., 2021, 121, e26339.
68. G. Bistoni, A. A. Auer and F. Neese, Chem. – Eur. J., 2017, 23, 865–873.
69. Q. Lu, F. Neese and G. Bistoni, Phys. Chem. Chem. Phys., 2019, 2019, 11569–11577.
70. L. P. Ho, A. Nasr, P. G. Jones, A. Altun, F. Neese, G. Bistoni and M. Tamm, Chem. – Eur. J., 2018, 24, 18922–18932.
71. R. Ghafarian Shirazi, F. Neese, D. A. Pantazis and G. Bistoni, J. Phys. Chem. A, 2019, 123, 5081–5090.
72. M. E. Beck, C. Riplinger, F. Neese and G. Bistoni, J. Comput. Chem., 2021, 42, 293–302.
73. S. Ghosh, F. Neese, R. Izsák and G. Bistoni, J. Chem. Theory Comput., 2021, 17, 3348–3359.
74. G. Bistoni, A. Altun, Z. Wang and F. Neese, Acc. Chem. Res., 2024, 57, 1411–1420.
75. A. Altun, I. F. Leach, F. Neese and G. Bistoni, ChemRxiv, 2024, preprint DOI:10.26434/chemrxiv-2024-fg6lj.
76. Q. Lu, F. Neese and G. Bistoni, Angew. Chem., Int. Ed., 2018, 57, 4760–4764.
77. Z. Wang and F. Neese, J. Chem. Phys., 2023, 158, 184102.
78. A. Altun, F. Neese and G. Bistoni, J. Chem. Theory Comput., 2019, 15, 5894–5907.
79. A. Altun, F. Neese and G. Bistoni, J. Chem. Theory Comput., 2022, 18, 2292–2307.
80. A. Altun, M. Garcia-Ratés, F. Neese and G. Bistoni, Chem. Sci., 2021, 12, 12785–12793.
81. J. Schümann, L. Ochmann, J. Becker, A. Altun, I. Harden, G. Bistoni and P. Schreiner, J. Am. Chem. Soc., 2023, 145, 2093–2097.
82. A. M. Preda, W. B. Schneider, D. Schaarschmidt, H. Lang, L. Mertens, A. A. Auer and M. Mehring, Dalton Trans., 2017, 46, 13492–13501.
83. C. Ganesamoorthy, S. Heimann, S. Hölscher, R. Haack, C. Wölper, G. Jansen and S. Schulz, Dalton Trans., 2017, 46, 9227–9234.
84. C. H. Borca, B. W. Bakr, L. A. Burns and C. D. Sherrill, J. Chem. Phys., 2019, 151, 144103.
85. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. De Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
86. P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. De Gironcoli, P. Delugas, R. A. Distasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H. Y. Ko, A. Kokalj, E. Kücükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H. V. Nguyen, A. Otero-De-La-Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu and S. Baroni, J. Phys. Condens. Matter, 2017, 29, 465901.
87. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
88. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
89. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
90. C. Zhang, J. Liu, Z. Ren, D. D. Koelling and B. N. Harmon, J. Phys. C: Solid State Phys., 1977, 10, 3107.
91. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188–5192.
92. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 73–78.
93. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1327.
94. F. Neese, F. Wennmohs, U. Becker and C. Riplinger, J. Chem. Phys., 2020, 152, 224108.
95. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2022, 12, e1606.
96. E. Rebolini, R. Izsák, S. S. Reine, T. Helgaker and T. B. Pedersen, J. Chem. Theory Comput., 2016, 12, 3514–3522.
97. S. Kossmann and F. Neese, Chem. Phys. Lett., 2009, 481, 240–243.
98. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
99. G. Bistoni, C. Riplinger, Y. Minenkov, L. Cavallo, A. A. Auer and F. Neese, J. Chem. Theory Comput., 2017, 13, 3220–3227.
100. S. F. Boys, Rev. Mod. Phys., 1960, 32, 296–299.
101. A. Altun, LEDAW – LED Analysis Wizard for Automating Local Energy Decomposition Analysis Using ORCA Outputs, 2024, https://github.com/ahmetaltunfatih/LEDAW.
102. A. Wuttke and R. A. Mata, J. Comput. Chem., 2017, 38, 15–23.

---

Footnote

† Electronic supplementary information (ESI) available: Detailed energetics and optimized coordinates of the structures (XLSX). See DOI: https://doi.org/10.1039/d4cp03697b

---