A combined density functional and coupled-cluster theory study on correlation-bound anions of perfluorinated compounds

Abstract

Perfluorinated cage molecules are hypothesized to have excellent electron-capture abilities, since they have the capacity to host an electron inside the carbon framework. The formed anions have been characterized as correlation-bound, meaning that they are unbound at the Hartree–Fock level. In this study, we assess the performance of density functional and coupled-cluster theory in the calculation of electron affinities of perfluorotetrahedrane, perfluorocubane, perfluoroadamantane and perfluorobenzene. We also characterize the anionic states using the electronic Fukui function and the electron localization function and investigate geometry changes upon electron attachment. To this end, the charge stabilization method is used for describing both the energy and spatial functions of metastable anionic states whenever they occur, and special attention is paid to how the electronic Fukui function and electron localization function differ for unbound anions compared to bound anions. For valence anions, density functional theory turns out to perform at a similar level as spin-scaled versions of second-order approximate coupled cluster singles and doubles theory, but it is not as accurate for describing correlation-bound anions. Only for the largest considered cage molecule, perfluoroadamantane, we found evidence for electron capture inside the cage, and this was only before the molecular geometry was allowed to relax. For anions that are described as metastable, the localization inside the cage became less pronounced.

---

1 Introduction

Many compounds form stable anions upon attachment of an excess electron. These can be categorized depending on the binding mechanism: they can be valence anions, dipole- or quadrupole-bound anions, or correlation-bound anions.^1,2 In this study, we focus on the latter type, where the excess electron is bound exclusively by correlation effects.^3–8 Therefore, Hartree–Fock theory, which lacks a description of electron correlation, will predict these anions to be unstable towards loss of an electron. A specific set of correlation-bound anions formed by perfluorinated cage molecules drew our attention because the excess electron is said to be bound inside the cage.^9–12 One of the earliest computational studies on these molecules was performed by Li and co-workers on perfluoroadamantane and its anion, focusing on equilibrium geometries, electron affinities and vibrational frequencies using four different density functional approximations.¹³ The same authors simultaneously performed a computational study on the perfluoroadamantyl radical and its anion, i.e. perfluoroadamantane with one fluorine atom removed (C₁₀F₁₅).¹⁴ The location of the excess electron was not studied, and presumably at that time, it was not yet postulated that the electron would reside inside the cage. Only later, in 2008, Irikura presented the hypothesis of σ-stellation: when the antibonding carbon–fluorine σ* orbitals overlap at a common point in space, they form a totally symmetric orbital that should be able to accommodate an electron.⁹ This was illustrated by a density functional theory study on ring and cage molecules with varying degrees of fluorination, with special attention to how the fluorination pattern influences the electron affinities of the molecules. The electron affinities were indeed found to increase with the number of overlapping σ* orbitals.⁹ Since then, interest in these molecules has increased and many authors have investigated them for their potential to trap electrons. Among these investigations are computational studies performed on a series of cage molecules, all employing density functional theory.^9,10,12,15 These comprise visualization of anionic singly occupied molecular orbitals (SOMOs)^9,15,16 or spin densities,^9,12 lowest unoccupied molecular orbitals of the neutral molecules¹² and/or Fukui functions,⁹ all suggesting that the excess electron would be mostly present inside the molecular cages. Recently, this even aroused interest in how nuclear spin–spin coupling is affected by the excess electron.^15,16 For some larger cages, such as perfluoroadamantane (C₁₀F₁₆) and perfluorododecahedrane (C₂₀F₂₀), electron localization functions (ELFs) showed a weak local maximum at the center of the anion, which had not been observed in the ELF for the neutral molecule.¹⁰ However, for the smaller cage molecules, this maximum was absent.^10,15 Moreover, integration of spin densities and Fukui functions over the cage interior yielded rather small electron populations, for both larger and smaller cage molecules, implying that the electron would be mostly outside the cage.^9,10

The synthesis of neutral perfluoroadamantane was first reported in the 1980s,^17,18 albeit not for the study of electron capture, the authors mentioned the use as a synthetic blood substitute as the rationale.¹⁸ To the best of our knowledge, synthesis of the perfluoroadamantane anion has not yet been reported. Recently, however, Sugiyama et al. succeeded in synthesizing perfluorocubane and its anion, driven by its hypothesized electron capture abilities.¹¹ They also reported evidence for the presence of the electron inside the cage, obtained by taking several steps: first, they measured the fluorine hyperfine coupling constants for the perfluorocubane anion they synthesized. Next, they calculated the corresponding values using density functional theory, finding that they were in excellent agreement with the experimental values. From this, they concluded that the calculated electronic state gives a good description of their synthesized anion. Finally, they assessed the spin density obtained from the calculations, which is mainly distributed evenly inside the cage, showing that the electron should be contained inside the cage. This gave rise to greater interest in this molecule and additional investigations. Computational studies on the perfluorocubane anion include the assessment of induced ring currents under an external magnetic field,^19,20 thermochemistry,²⁰ and topological analysis of the electron density.²⁰ The corresponding results all indicate that perfluorocubane is a good electron acceptor.^19,20 Recently, Gu and Jin even investigated the possibility of having two electrons trapped in the perfluorocubane cage by using density functional theory to calculate the molecular orbitals of a complex formed by coating perfluorocubane with one or two metal atoms. The metal atoms were found to each transfer an electron to the interior of the cage.²¹ We would also like to mention here the recent review of Krafft and Riess,²² presenting an extensive overview of what has been done and what is known in the field of cage molecule research, both experimentally and computationally. They also expressed criticism by noting that if an electron can readily enter the cage, it can also easily exit it again. Furthermore, they emphasized that experimental evidence has been limited so far. This leads to one objective of the current study, which is to investigate these anions using a broad range of quantum chemical methods, with special attention for their correlation-bound character and the location of the excess electron. Recently, the current authors computed Fukui functions and electron localization functions of metastable anions²³ by applying the charge stabilization and extrapolation method.²⁴ This approach will also be applied here to metastable anions whenever they occur. In particular, we are curious whether the results for the metastable anions predicted using Hartree–Fock theory are in line with those for the bound anions predicted using density functional and coupled-cluster theory.

For the current study, we selected three cage molecules: perfluorotetrahedrane (C₄F₄), perfluorocubane (C₈F₈), and perfluoroadamantane (C₁₀F₁₆). First, we assess the performance of different methods for the calculation of electron affinities of these molecules, where we also include perfluorobenzene (C₆F₆). Although this molecule is not a cage, it forms a correlation-bound anion. Additionally, it shows significant geometry changes upon electron attachment, while turning into a valence anion. The equilibrium geometry of the neutral molecule is planar and displays D_6h symmetry. After electron attachment and geometry relaxation, the resulting molecule has a buckled structure with C_2v symmetry.^6,25–27 We will therefore pay particular attention to how different methods perform in the description of the anionic structures and states of this molecule, taking into account that some methods predict an unbound anion while others do not. Reports on the performance of second-order approximate coupled-cluster methods²⁸ and density functional theory methods²⁹ for correlation-bound anions are limited, providing motivation for the method assessment presented here. For the cage molecules C₄F₄, C₈F₈ and C₁₀F₁₆, we characterize the electronic structures of the anions of the cage molecules, using the electronic Fukui functions and the electron localization function, with special attention to where the electron is localized. Additionally, we monitor the performance of HF theory compared to the density functional methods, as well as whether densities and derived quantities for unbound states differ significantly from those for bound states. Whenever unbound anions occur, they are described using the extension to the charge stabilization and extrapolation method²⁴ recently proposed by the current authors.^23,30 It had previously been applied to the temporary anions of chlorinated ethene derivatives that are subject to dissociative electron attachment.^23,30 The molecules studied here provide a different set of systems with different characteristics, and hence they are ideal for further assessment of our method. The details of the performed calculations can be found in Section 2. In Section 3, the results for electron affinities are presented together with a comparison of the used methods and basis sets. Section 4 discusses the electronic structures of the studied anions in more detail. Finally, we present the conclusions of the study in Section 5.

2 Computational details

2.1 Electronic structure calculations

The geometries of neutral and anionic species were optimized at the ωB97M-V³¹/aug-cc-pVTZ^32,33 level of theory. The optimization tolerances for the maximum gradient component, the maximum atomic displacement, and the maximum absolute energy change were set to 10⁻⁶, 10⁻⁶ and 10⁻⁸, respectively. The SCF convergence criterion was set to 10⁻⁸. The coordinates of the optimized geometries are available in the SI. For neutral C₄F₄, no stable geometry was found.

Single-point energy calculations were performed using Hartree–Fock (HF) theory,^34–36 three different density functionals (PBE,³⁷ M06³⁸ and ωB97M-V³¹), the second-order approximate coupled-cluster singles and doubles theory with the resolution-of-the-identity (RI-CC2)³⁹ and two spin-scaled versions thereof (SOS-RI-CC2⁴⁰ and SCS-RI-CC2⁴¹), and finally equation-of-motion coupled-cluster theory with single and double excitations (EOM-CCSD).⁴² The employed density functional approximations belong to the categories of generalized gradient approximation (GGA), global hybrid and range-separated hybrid, respectively. PBE was selected since it is a popular choice of GGA⁴³ and M06 because it was found to perform well for the electron affinities of metastable anions.⁴⁴ ωB97M-V was chosen for its good performance as a “workhorse functional”⁴⁵ and because the use of functionals with exact long-range exchange was proposed as a solution to underestimation of the binding energy of loosely bound electrons.^46,47 For DFT numerical quadrature, the SG-2 grid⁴⁸ was applied.

The two used basis sets are aug-cc-pVTZ and aug-cc-pVTZ with 7 additional diffuse s and p shells. This latter basis set is inspired by the ones used by Voora and Jordan for the description of the anion of perfluorobenzene⁶ and by Paran and co-workers for assessing the performance of CC2 methods for describing non-valence anions.²⁸ The two basis sets differ from one another by the position of the additional shells: in the study by Voora and Jordan, they were positioned only at the center of the molecule, while in that by Paran and co-workers and in our study, they were positioned at every atomic nucleus instead to ensure that the basis set is nearly complete. For every new shell, the exponent was determined to be half of the exponent of the previously most diffuse s or p shell. As is typical when many diffuse basis shells are used, linear dependencies can occur. These are detected and removed before the SCF equations are solved.

In the case of metastable anions, the charge stabilization method²⁴ was applied according to the procedure explained in ref. 30. This method comprises scaling the nuclear charges of the molecule to create an additional attractive potential for stabilizing the metastable state. All electronic energy values corresponding to bound states are subsequently extrapolated to obtain the electronic energy of the unperturbed metastable state.²⁴ While methods exist that yield a complex-valued electronic energy for such states, whose imaginary part describes the rate of decay,⁴⁹ the charge stabilization method employed here yields real-valued energies.

The adiabatic electron affinity (AEA) without zero-point energy (ZPE) correction was taken as the electronic energy difference between the neutral and anionic states at the respective equilibrium geometries,

AEA = E(neutral at neutral eq. geom.) − E(anion at anion eq. geom.) (1)

The vertical electron affinity (VEA) of the neutral compound was calculated as the electronic energy difference between the neutral and anionic states in the equilibrium geometry of the neutral state:

VEA = E(neutral at neutral eq. geom.) − E(anion at neutral eq. geom.) (2)

Finally, the vertical detachment energy (VDE) of the anionic compound was determined as the electronic energy difference between the two states in the equilibrium geometry of the anionic state:

VDE = E(neutral at anion eq. geom.) − E(anion at anion eq. geom.) (3)

For the HF and DFT calculations, VEA and VDE were obtained via two separate calculations using the same method and basis set. For the EOM-CCSD and RI-CC2 calculations, they were obtained in a single calculation, using the electron attachment variant (EOM-EA)⁵⁰ with the ground states of the neutral molecules as CC reference states.

Electronic structure calculations were performed using the Q-Chem 6.0.2 software package.⁵¹ Energies reported for C₈F₈ at the EOM-CCSD/aug-cc-pVTZ level of theory were obtained using the CFOUR 2.1 software package.⁵²

2.2 Spatial functions

For metastable anions, spatial functions were obtained from charge stabilization and extrapolation, by applying our recently developed procedure.²³ Through theoretical prediction of dissociative electron attachment to chlorinated ethene derivatives, we showed that the obtained extrapolated spatial functions are reliable,²³ justifying the use of the same procedure in the current research. This allows for direct comparison between spatial functions for metastable anions and bound anions.

Formatted checkpoint files were produced using Q-Chem from the HF and DFT single-point energy calculations. These were converted to extended wavefunction (wfx) files using the conversion tool included in the AIMAll package.⁵³ From these wfx files, highest occupied molecular orbitals (HOMOs), electron densities and electron localization functions (ELFs)^54,55 were obtained using the Multiwfn package⁵⁶ and stored as Gaussian-type cube files. Electronic Fukui functions f⁺(r)⁵⁷ were calculated by subtracting the density of the neutral molecule (describing N₀ electrons) from that of the anion (describing N₀ + 1 electrons):

f⁺(r) = ρ_N0+1(r) − ρ_N0(r) (4)

Using the property that the integral of the electron density over the complete space yields the number of electrons it represents, it is derived that the integral of the Fukui function equals 1. Subsequently, by integrating the Fukui function only over the cage interior, the resulting value can be interpreted as the fraction of the excess electron residing inside the cage. These integrations were performed numerically using the values stored in the produced cube files. In practice, due to small numerical inaccuracies, the value of the integral over the complete file is not exactly equal to 1 but varies between 0.96 and 1.03. Therefore, the value obtained from integrating over the cage interior was scaled accordingly by dividing it by the integral over the complete file.

To ensure that calculations performed using different methods and/or different basis sets result in the same electronic states, the anionic HOMOs obtained from HF theory and DFT were compared to Dyson orbitals calculated from the EOM-CCSD and CCSD wave functions.^58,59 These Dyson orbitals describe the overlap between the N₀ and N₀ + 1 wavefunctions:

The radial distribution functions shown in Fig. 1 were calculated using the qcip_tools package⁶⁰ for Python.

Fig. 1 Left panel: highest occupied molecular orbital of the anionic state of C₆F₆ and the Dyson orbital for electron attachment to neutral C₆F₆ (isovalue = 0.001 a.u.), at equilibrium geometry of the neutral compound. Obtained using M06/aug-cc-pVTZ+7s7p (by charge stabilization and extrapolation), ωB97M-V (by charge stabilization and extrapolation), and EOM-EA-CCSD/aug-cc-pVTZ+7s7p. Right panel: radial distribution functions of the visualized orbitals.

3 Attachment and detachment energies

The obtained VEAs, AEAs (without ZPE correction) and VDEs are reported in Table 1. First, the influence of adding more diffuse shells to the basis set is discussed. The degree to which this alters the obtained VEA, AEA and VDE values appears to depend on the described molecule and property. In general, the changes in AEA values are smaller than those in VEA and VDE values. Possibly, this results from a cancellation of errors, since for the calculation of the AEA, two geometry optimizations are required (one for the neutral compound and the other for the anion), while for the VEA and the VDE, only one optimization is needed. The smallest energy changes upon increasing the basis set are found for AEA and VDE values for C₁₀F₁₆, regardless of which method is used for their calculation. In contrast, the largest energy changes are observed for VEA values for C₆F₆ and VDE values for C₄F₄, again independent of the used method. In many cases, these changes are larger than the conventional threshold of chemical accuracy, i.e. 1 kcal mol⁻¹, corresponding to 0.0434 eV, which is a sign that the results are not yet converged with respect to the basis set. Furthermore, in some cases, the character of the anion changes from unbound to bound as the basis set is increased; this can be seen for the VEA for C₁₀F₁₆ calculated using SOS-RI-CC2, the VEA for C₆F₆ calculated using SCS-RI-CC2, SOS-RI-CC2 and EOM-CCSD, and for the VDE for C₄F₄ calculated using M06. This shows that additional shells are essential for binding the excess electron. Therefore, it is concluded that the predefined standard basis set (aug-cc-pVTZ) is insufficient for accurately describing the attachment and detachment energies of the studied compounds. This is in agreement with previous research based on CC methods,^28,61 and our work now shows that DFT methods have similar requirements regarding the basis set.

Table 1 Attachment and detachment energies (eV) computed using different methods specified in the first column (aug = aug-cc-pVTZ, aug+7 = aug-cc-pVTZ+7s7p). Values in italics indicate unbound anions. Bold indicates our best values for the attachment and detachment energies of the studied compounds

Method\basis set	C4F4		C8F8		C10F16		C6F6
	aug	aug+7	aug	aug+7	aug	aug+7	aug	aug+7
a No values due to the absence of equilibrium geometry for neutral C4F4. b No values due to SCF convergence issues with PBE/aug-cc-pVTZ+7s7p. c No values due to unphysical electron densities obtained using M06/aug-cc-pVTZ+7s7p. d No values due to too high memory requirements. e Unbound, no values because no charge stabilization was applied.
Vertical electron affinities
HF			−0.4349	−0.4058	−1.4170	−1.4291	−1.0797	−0.3305
PBE			1.3168		0.6998		0.1687
M06			0.8666	0.9414	0.3899		−0.4636	−0.4710
ωB97M-V			0.8516	0.8623	0.0606	0.1074	−0.3822	−0.2037
RI-CC2			1.6820	1.6861	0.9302	0.9384	0.2757	0.3785
SCS-RI-CC2			1.0789	1.0871	0.2541	0.2905	Unbound	0.0855
SOS-RI-CC2			0.7868	0.8004	Unbound	0.0473	Unbound	0.0147
EOM-CCSD			1.0280				Unbound	0.0871
Adiabatic electron affinities (without ZPE)
HF			−0.4987	−0.4876	−0.2577	−0.2567	−0.7080	−0.7071
PBE			0.0648		1.2687	1.2700	0.6907
M06			1.0633	1.0678	0.8463		0.2628	0.2662
ωB97M-V			1.1644	1.1672	0.8243	0.8249	0.2973	0.3028
RI-CC2			2.1185	2.1204	1.9853	1.9878	1.0421	1.0458
SCS-RI-CC2			1.4657	1.4686	1.2407	1.2430	0.6043	0.6101
SOS-RI-CC2			1.1437	1.1477	0.8693	0.8718	0.3947	0.4021
EOM-CCSD			1.2427				0.2889	0.2963
Vertical detachment energies
HF	−0.7869	−0.4684	0.2281	0.2393	3.2382	3.2385	0.8400	0.8423
PBE	0.4390		1.9986		3.2890	3.2896	1.5493
M06	−0.0179		1.6009	1.6061	3.5470	3.5465	1.4841	1.4886
ωB97M-V	0.0819	0.1686	1.5250	1.5282	3.4025	3.4026	1.3713	1.3778
RI-CC2	0.5889	0.6230	2.4077	2.4099	4.1379	4.1389	2.0096	2.0130
SCS-RI-CC2	0.1893	0.2678	1.7609	1.7640	3.4049	3.4057	1.5675	1.5729
SOS-RI-CC2	0.0071	0.1298	1.4402	1.4446	3.0344	3.0354	1.3543	1.3612
EOM-CCSD	0.2111	0.2902	1.6896				1.4175	1.4250

---

Next, the performance of the different employed methods is considered. Wherever possible, EOM-CCSD is used as the reference method. Due to the small number of data points, a rigorous analysis by means of, e.g., root-mean-square deviation is not possible. Instead, the analysis is based on observations made from Table 1. A first observation is that HF theory often describes anions as unbound where other methods describe them as bound, as seen from the VEA and AEA values for C₈F₈, C₁₀F₁₆ and C₆F₆ as well as the VDE values for C₄F₄. This is an indication of the correlation-bound character of the anions. However, as HF theory does predict the VDE values to be positive, these anions are no longer correlation-bound after structural relaxation. We would like to note here that this enables us to obtain the AEA values without the use of the charge stabilization method, even when these are negative; this is possible since the anion is a stable species at its equilibrium geometry.

The usefulness of SCS- and SOS-RI-CC2 is highlighted by the fact that these methods produce results for cases where EOM-CCSD becomes unfeasible due to large system size (indicated by footnote “d” in Table 1). In general, the values obtained from the spin-scaled RI-CC2 methods (SCS-RI-CC2 and SOS-RI-CC2) are closer to the EOM-CCSD value than those obtained from the (unscaled) RI-CC2 method. Our findings regarding the accuracy of VEA values obtained using spin-scaled RI-CC2 methods are in line with what was reported in an earlier study on non-valence anions: they always provide an improvement over unscaled RI-CC2 methods, even though in some cases, the final result is not accurate, while for others, it is in excellent agreement with EOM-CCSD.²⁸ Notably, SOS-RI-CC2 performs better for the calculation of AEAs, while for vertical quantities (VEA and VDE), SCS-RI-CC2 performs best. The cause of these trends remains unclear, but we note that our results for VEA match the trends seen in earlier work on valence, dipole-bound and correlation-bound anions.²⁸ Yet, the opposite trend, with SOS outperforming SCS, was observed for the closely related ADC(2) method applied to VEAs of valence anions.⁶² Further data allowing us to compare the performance of SCS and SOS for computing adiabatic quantities are, to our best knowledge, not yet available. More cases need to be studied to see whether this trend extends beyond the current research.

Our best values for the attachment and detachment energies are determined as follows: results from EOM-CCSD are used as much as possible, with preference for the more diffuse basis set. In those cases where no EOM-CCSD results are available, SCS-RI-CC2 (for VEA and VDE) or SOS-RI-CC2 (for the AEA) in combination with the more diffuse basis set is preferred. The corresponding values are shown in Table 1 in bold.

Concerning the density functional approximations, the results imply that the hybrid density functional approximations M06 (global hybrid) and ωB97M-V (range-separated hybrid) perform better than the generalized gradient approximation PBE. This latter functional also suffered from convergence issues when used in combination with the aug-cc-pVTZ+7s7p basis set for describing the anionic states. Furthermore, the results from M06 and ωB97M-V are often very similar to those from SCS- and SOS-RI-CC2, in particular for AEA and VDE values. For calculating these properties, the energy of the anionic state at its equilibrium geometry is required, and this is precisely when the anion has valence character rather than correlation-bound character. This indicates that DFT methods exhibit poor performance in describing correlation-bound anions, while for valence anions, in particular using range-separated hybrid functionals, it is possible to obtain accurate results at lower cost as compared to CC2 methods. This is in line with earlier reports stating that range-separated functionals perform well for the description of anions,^29,30,63 while they fail in describing the electron binding energy of the correlation-bound anion of C₆F₆.²⁹ Furthermore, functionals with exact exchange in the long range were proposed as a solution for avoiding underestimation of the electron binding energies as the basis set increases,⁴⁶ and our findings confirm this. Additionally, M06 suffered from producing unphysical electron densities when combined with the more diffuse basis set. Such densities were observed for anionic C₄F₄ at the equilibrium geometry of the anion and anionic C₁₀F₁₆ at the equilibrium geometry of the neutral molecule. For both of these, the electron density did not integrate to the number of electrons but instead lacked about 0.5 electron. M06 is known to be more sensitive to the quality of the integration grid, posing problems in describing anions using diffuse basis sets. For the case of the anions described here, it is thus found that ωB97M-V performs better than M06.

Finally, the results for C₆F₆ are discussed in more detail. Our reference value of 0.087 eV, obtained at the EOM-CCSD/aug-cc-pVTZ+7s7sp level of theory, is close to the value of 0.135 eV reported by Voora and Jordan obtained using a similar method and basis set.⁶ The discrepancy between the two values is in line with the trend of EOM-CCSD yielding lower VDE values than EOM-CCSD(2).⁶¹ When the more compact basis set is employed, it is seen that nearly all methods except PBE and RI-CC2 describe the anion as unbound. This is related to what was mentioned above, which is that additional diffuse shells are needed for binding the excess electron. Yet, after expanding the basis set, negative values are still found using Hartree–Fock theory, M06 and ωB97M-V. For the first method, this is expected, but from the last two methods, it would be expected that a bound anion is found. Comparison of the extrapolated M06 and ωB97M-V anionic HOMOs with the EOM-EA-CCSD/aug-cc-pVTZ+7s7p Dyson orbital for electron attachment, both visually and via the radial distribution function, reveals that both HOMOs are similar to the EOM-CCSD Dyson orbital at a qualitative level (Fig. 1). However, the Dyson orbital is more diffuse. We thus conclude that M06 and ωB97M-V, while being able to detect the electronic state for the correlation-bound anion of C₆F₆, do not describe the orbital hosting the excess electron accurately enough to make the anion bound. Given the fact that describing this anion is so challenging, it could be worthwhile to include it in benchmarks for new density functional approximations.

4 Characterization of the electronic states and molecular geometries

We performed a more in-depth characterization of the electronic states of the cage molecules C₄F₄, C₈F₈ and C₁₀F₁₆. Visualization of the LUMOs of the neutral compounds (see Fig. S2 in the SI) reveals that they are composed of the overlap of the carbon–fluorine σ* orbitals, in line with the σ-stellation hypothesis.⁹ This allows the electron to reside inside the cage, which we assess in more detail by analyzing the computed electronic Fukui functions and electron localization functions. It is important to note that these belong to the same irreducible representation as the orbital to which the electron is attached. From the Fukui functions, we computed which fraction of the excess electron resides in the cage interior. The obtained percentages are plotted in Fig. 2. The electron localization functions were used to examine whether a maximum of the electron localization is detected at the cage center. In the following subsections, the results are discussed for each considered molecule separately.

Fig. 2 Fraction of excess electron inside the cage, determined by integration of the electronic Fukui function. For details, see the text. Stars indicate results for vertically unbound anions (aug = aug-cc-pVTZ and aug+7 = aug-cc-pVTZ+7s7p).

4.1 Perfluorotetrahedrane

C₄F₄, with the trivial name “perfluorotetrahedrane”, is a hypothetical molecule with a tetrahedral structure belonging to the T_d point group. There is no consensus about the equilibrium geometry of the neutral molecule. Some authors succeeded in finding a tetrahedral equilibrium geometry using the semi-empirical PM3 method,⁶⁴ HF theory,⁹ the B3LYP* density functional approximation that was developed for the description of transition metals¹² or the conventional B3LYP density functional approximation.¹⁵ However, other authors report that the tetrahedral geometry does not represent a minimum and shows vibrations with imaginary frequencies. Among these are optimizations performed using the B3LYP density functional^10,65,66 or MP2.⁶⁵ Regardless of whether or not the neutral molecule displayed a tetrahedral equilibrium structure, in all of these studies, the anion was found to have a T_d equilibrium geometry.^9,10,12,15 In the current study, our attempt to optimize the neutral state at the ωB97M-V/aug-cc-pVTZ level of theory led to a four-membered ring structure. The T_d geometry of minimal energy also turned out to be a transition state rather than an equilibrium geometry. For the anionic state, we did obtain a T_d structure as equilibrium geometry, where the electronic state is of A₁ symmetry. An investigation of the underlying mechanism for the instability of the neutral molecule is reported in the SI, where it is shown that the excess electron contributes to holding the molecule together.

The obtained electron fractions from integration over the cage are in line with the observations from the visualized Fukui function (Fig. S1 and S2 in the SI): they are not larger than 3.35%, meaning that the excess electron is situated at least 96% outside of the cage. Note that the fraction for HF/aug-cc-pVTZ+7s7p is even slightly negative (−0.27%), which is interpreted as (albeit very limited) electron depletion from the cage interior upon electron attachment.

The anionic ELFs did not show a maximum at the center of the cage (see Fig. 3), in line with earlier findings by Berski et al.¹⁰ Rather, they show a minimum, regardless of whether the anion is described as bound or unbound. Additionally, the difference between anionic and neutral ELFs showed that there is a very small increase in electron localization of at most 0.11 a.u. inside the cage, for all combinations of methods and basis sets except for HF/aug-cc-pVTZ+7s7p. The absence of this increase for this latter case is logical, given the electron depletion already observed from the integrated Fukui functions. For all combinations of methods and basis sets, there is also increased electron localization outside of the cage, mostly of similar magnitude to the increase at the center of the cage. Specifically, for HF theory, however, the increase outside of the cage is larger (up to 0.21 a.u.) and even leads to an amount of localization that is visible in the isocontour plots (indicated by arrows in Fig. 3). These are cases where the anion is described as unbound and the results were obtained using charge stabilization and extrapolation. Nevertheless, such a substantial amount of localization outside of the cage is not observed for M06/aug-cc-pVTZ, while it also describes the anion as unbound. It is thus most likely a consequence of a poor description by HF theory, rather than the use of the charge stabilization method. Combining all findings, it is concluded that there is no evidence for the presence of the electron inside the cage. Furthermore, there is no indication that the use of charge stabilization introduces deviations in the results. The observed differences between methods and the fact that certain methods describe the states as unbound where others describe them as bound seem to originate from the used level of theory itself.

Fig. 3 Isocontour plots of the anionic ELF at the equilibrium geometry of the C₄F₄ anion for different methods and basis sets (specified in the figure). The ELFs were plotted in a plane through the center of the molecule and containing two C–F bonds (see the inset at bottom left). The stars indicate the results for vertically unbound anions that were obtained using charge stabilization. The arrows indicate localization observed outside the cage.

4.2 Perfluorocubane

The structure of C₈F₈, with the trivial name “perfluorocubane”, is expected to be cubic in analogy to the cubic molecule cubane (C₈H₈).⁶⁷ This is the case for the equilibrium structure of the neutral molecule: the carbon atoms are arranged according to the corners of a cube, and the molecule belongs to the O_h point group. After electron attachment and structural relaxation, the new structure is still cubic but with slightly shortened carbon–carbon bonds and elongated carbon–fluorine bonds. These cubic structures were also reported by other authors before, both experimentally¹¹ and computationally.^{9,10,12,15,16,19} The electron attaches to an orbital of a_1g symmetry, explaining the absence of a Jahn–Teller effect. A pseudo-Jahn–Teller effect could still occur, but is not observed here.

Sugiyama and co-workers reported evidence for the presence of the electron inside the cage, based on the computed spin density.¹¹ However, from other theoretical studies, the electron fraction inside the cage is reported to be rather small.^9,10 This is also the case in the current study: the obtained fractions are all below 6.5%, and here as well, the slightly lower values for the more diffuse basis set indicate that the excess electron can occupy a broader region in space. Furthermore, the values do not change significantly upon geometry relaxation except for the cases described by HF theory. Before and after structural relaxation, the values obtained from HF/aug-cc-pVTZ are equal to 5.14 and 6.46%, respectively, while those obtained from HF/aug-cc-pVTZ+7s7p are equal to 4.39 and 6.24%, respectively. Furthermore, no maximum of the anionic ELFs is observed in the center of the cage, for any combination of method, basis set and geometry, confirming earlier findings.¹⁰ Instead, local minima are detected (see Fig. 4). A closer look at the Fukui function (see Fig. S3 in the SI) and the ΔELF between neutral and anionic states reveals that the excess electron is mostly located outside of the cage, near the fluorine atoms. Furthermore, upon adding diffuse shells, the obtained fractions slightly decrease. This is interpreted as that the more diffuse basis set allows the electron to occupy a more extended region in space and be less confined to the center of the molecule. In analogy to the study by Sugiyama et al.,¹¹ we also calculated the spin density for the anionic state at the equilibrium structures of both the neutral and anionic states (see Fig. S3 in the SI). The spin density is indeed located mainly inside the cage, but it is also located in the vicinity of the fluorine atoms. Hence, we conclude that the spin density alone is not sufficient to determine the position of the excess electron.

Fig. 4 Isocontour plots of the anionic ELF at the equilibrium geometry of the C₈F₈ anion, for different methods and basis sets (specified in the figure). The ELFs were plotted in a plane through the center of the molecule and containing two opposite C–C bonds and four C–F bonds (see the inset at bottom left).

4.3 Perfluoroadamantane

Perfluoroadamantane, C₁₀F₁₆, is the largest molecule included in our study set. The equilibrium geometry of the neutral molecule belongs to the T_d point group. In our current study, we found a C_3v structure with one elongated carbon–fluorine bond (1.98 Å versus 1.35 Å) as geometry of minimal energy for the anion. A further discussion on these structures is included in the SI.

Anionic perfluoroadamantane at its equilibrium geometry yields a larger fraction of the excess electron contained in the cage interior, as compared to the other molecules (see Fig. 2). At the equilibrium geometry of the neutral compound, the fractions vary between 3.16 and 5.62%, while at that of the anion, they lie between 6.48 and 8.66%. These latter values are in line with what Irikura obtained for the integration of the difference between anionic and neutral electron densities (6.4%).⁹ Yet, this still indicates that there is a higher fraction of the electron outside the cage than inside.

Looking at the anionic ELFs obtained at the equilibrium geometry of the neutral compound (see Fig. 5), a maximum is observed in the cage center, albeit very weak. The maximum has a value between 0.10 and 0.17 a.u., depending on the used method and basis set. The largest value of 0.17 a.u. is very close to the value of 0.176 a.u. reported by Berski and co-workers.¹⁰ For reference, other maxima such as those found in bonds or lone pairs take values between 0.7 and 0.9 a.u. We also determined the distance between the ELF maximum and the nearest other ELF critical points, as well as the corresponding values. The nearest critical points are all saddle points, and they are located less than 0.5 Å away from the central maximum. The ELF values at these saddle points are also nearly the same as those of the maximum, demonstrating that the maxima at the cage center are not very pronounced. The exact ELF values at maxima and saddle points and the distances between them are listed in the SI. Additionally, a maximum at the cage center was observed only when the anion was described as stable, i.e. when a DFT method was employed. When HF theory is used, the anion appears as metastable and the maximum in the anionic ELF is absent, even though the values of the integrated Fukui function are similar to those for DFT methods. We interpret this as the electron still being partly inside the cage, but less localized in the center. To determine whether this could be an artifact of the extrapolation procedure, we computed an extrapolated ELF using ωB97M-V/aug-cc-pVTZ, for which the anion is in fact already bound without the use of charge stabilization. This means that we scale the nuclear charges and extrapolate quantities of interest in the same way as for an unbound state, with the difference that we can now compare to the quantity computed directly for the unperturbed bound system. The extrapolated ELF is nearly identical to the directly computed ELF, with the absolute difference between them not exceeding 0.006 a.u. Both show an ELF maximum at the same location in the cage center, with values of 0.150 a.u. (extrapolated ELF) and 0.144 a.u. (directly computed ELF). This indicates that the absence of a maximum for the states described by HF theory stems from the use of HF theory itself rather than the use of the charge stabilization and extrapolation procedure. This is another indication that HF theory is not able to accurately describe these states.

Fig. 5 Isocontour plots of the anionic ELF at the equilibrium geometry of neutral C₁₀F₁₆, for different methods and basis sets (specified in the figure). The ELFs were plotted in a plane through the center of the molecule and containing four carbon atoms (see the inset at bottom left). Stars indicate results for vertically unbound anions obtained using charge stabilization, and the arrows indicate the very weak localization in the cage center.

At the equilibrium geometry of the anion, the maximum in the anionic ELFs disappears (see Fig. S4 in the SI). Instead, increased localization is found near the elongated carbon–fluorine bond, akin to localization in an antibonding σ* orbital. This is best seen when the ΔELF, the difference in the ELF occurring on electron attachment, is plotted next to the HOMO of the anion (see Fig. 6). The HOMO is still situated partly inside the cage, most likely explaining the relatively large electron fraction value obtained for this system. We also calculated the ELF for the structure at which the anion has the lowest energy within the T_d point group. The increased localization akin to a σ* orbital is no longer present, but the maximum reappears, although again only when the anion is bound. Apparently, the tetrahedral symmetry prevents the excess electron to be localized in one single bond, as this would imply symmetry breaking.

Fig. 6 Left and middle panels: isosurface plots of ΔELF for C₁₀F₁₆ at the equilibrium geometries of the neutral (left) and anionic (middle) states (isovalue = +0.05 a.u.). Right panel: isosurface plot of the highest occupied molecular orbital of C₁₀F₁₆⁻ at the equilibrium geometry of the anionic state (isovalue = 0.03 a.u.). The elongated carbon–fluorine bond is shown in the bottom left corner, which is computed at the ωB97M-V/aug-cc-pVTZ+7s7p level of theory.

We summarize our findings for C₁₀F₁₆ as follows: upon electron attachment, initially, a highly symmetric anion is formed which contains the electron partly inside the cage. This also occurs when the anion is described as unbound, but in these cases, the excess electron is less localized at the cage center. According to our geometry optimizations using ωB97M-V/aug-cc-pVTZ, the anionic potential energy surface does not have a minimum in the tetrahedral structure and relaxes to a structure of C_3v symmetry, where one carbon–fluorine bond shows substantial elongation. This is accompanied by a change in the location of the excess electron, which after relaxation is mainly situated at the elongated bond.

5 Conclusions

We investigated four perfluorinated molecules and their correlation-bound anions. The aims were two-fold: assessing different methods and basis sets for the description of attachment and detachments energies for these compounds, and gaining insight into the localization of the excess electron in three different cage molecules and how this is influenced by whether the anion is bound. This work is one of the first applications of the charge stabilization method to spatial functions and highlights how electron densities and electron localization functions for metastable states can now be treated on an equal footing with those for bound states.

For the RI-CC2 methods, there are indications that spin-scaling improves the accuracy of the obtained electron affinities and detachment energies. Additionally, two of the density functional approximations, M06 and ωB97M-V, were found to exhibit similar performance to the spin-scaled RI-CC2 methods for calculating the attachment and detachment energies of valence anions, meaning that similar accuracy can be achieved using less computational resources. Confirming previous research,²⁹ obtaining an accurate description of correlation-bound anions remains hard, even though employing the charge stabilization method enables M06 and ωB97M-V calculations to converge towards a state qualitatively corresponding to the correlation-bound anion. Furthermore, very diffuse shells in the basis set turned out to be essential for accurate calculation of the attachment and detachment energies. Yet, from the electronic Fukui functions and the electron localization functions, it became clear that these additional shells are not necessary for obtaining a qualitative description of the electronic state itself. Moreover, the use of the charge stabilization method turned out not to introduce inaccuracies; rather, observed inaccuracies were due to the use of HF theory, whose performance remained poor for these properties as well. These findings all highlight the difficulties associated with studying correlation-bound anions.

For the cage molecules, Irikura's hypothesis of overlapping carbon–fluorine σ* orbitals (σ-stellation)⁹ was assessed. The LUMOs of the neutral compounds, as well as the observed changes in bond length following electron attachment to perfluorocubane (C₈F₈), are in agreement with this hypothesis. Yet, for the two smallest cage molecules, perfluorotetrahedrane (C₄F₄) and perfluorocubane, we found no evidence for the presence of the excess electron inside the cage. For the largest molecule, perfluoroadamantane (C₁₀F₁₆), we found that initially, the electron can be localized in the cage except when it is described as unbound. Upon geometry relaxation, one of the carbon–fluorine bonds becomes elongated and the electron is found to be mainly localized in regions akin to an antibonding σ* orbital of this bond. Following this trend, anions of larger cage molecules such as perfluorododecahedrane (C₂₀F₂₀) would thus be expected to have an even more pronounced tendency for hosting an electron inside the cage.^10,12 Future work could thus include applying our analyses to these molecules as well.

Conflicts of interest

The authors have no conflicts to disclose.

Data availability

The data supporting the findings of this study have been included as part of the supplementary information (SI). They contains the coordinates of the optimized equilibrium structures of all compounds, and corresponding discussions for C₄F₄ and C₁₀F₁₆, an overview of the data points included for the extrapolation procedures, selected figures of the lowest unoccupied molecular orbitals of the neutral cage molecules, the numerical values for the integral of the Fukui functions over the cage interiors, selected figures of the Fukui function and the spin density for C₈F₈, the numerical values of the ELF maximum in the anionic ELF of C₁₀F₁₆ and isocontour plots of this ELF. See DOI: https://doi.org/10.1039/d5cp00960j.

Acknowledgements

C. T. thanks the Research Foundation-Flanders (FWO) for the received PhD fellowship (grant no. 1103924N). T.-C. J. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant/agreement no. 851766), and the KU Leuven internal funds (grant C14/22/083). F. D. P. acknowledges the Strategic Research Program funding from the VUB. All authors thank the anonymous reviewers for their valuable comments that led to significant improvement of the manuscript. Fig. 1 and 3–6 as well as the figures in the SI were created using VMD (for LinuxAMD64, version 1.9.4a51)⁶⁸ and GaussView 6.⁶⁹

References

1. J. Simons, J. Phys. Chem. A, 2008, 112(29), 6401–6511.
2. J. Simons, J. Phys. Chem. A, 2023, 127(18), 3940–3957.
3. T. Sommerfeld, B. Bhattarai, V. P. Vysotskiy and L. S. Cederbaum, J. Chem. Phys., 2010, 133(11), 114301.
4. V. G. Bezchastnov, V. P. Vysotskiy and L. S. Cederbaum, Phys. Rev. Lett., 2011, 107, 133401.
5. V. K. Voora, L. S. Cederbaum and K. D. Jordan, J. Phys. Chem. Lett., 2013, 4(6), 849–853.
6. V. K. Voora and K. D. Jordan, J. Phys. Chem. A, 2014, 118(35), 7201–7205.
7. V. K. Voora and K. D. Jordan, J. Phys. Chem. Lett., 2015, 6(20), 3994–3997.
8. V. K. Voora, A. Kairalapova, T. Sommerfeld and K. D. Jordan, J. Chem. Phys., 2017, 147(21), 214114.
9. K. K. Irikura, J. Phys. Chem. A, 2008, 112(5), 983–988.
10. S. Berski, A. J. Gordon and Z. Latajka, J. Phys. Chem. A, 2014, 118(23), 4147–4156.
11. M. Sugiyama, M. Akiyama, Y. Yonezawa, K. Komaguchi, M. Higashi, K. Nozaki and T. Okazoe, Science, 2022, 377(6607), 756–759.
12. A. Ghosh and J. Conradie, ACS Omega, 2023, 8(5), 4972–4975.
13. Q.-S. Li, X.-J. Feng, Y. Xie and H. F. Schaefer, J. Phys. Chem. A, 2005, 109(7), 1454–1457.
14. X.-J. Feng, Q.-S. Li, Y. Xie and H. F. Schaefer, J. Chem. Theory Comput., 2005, 1(2), 279–285.
15. X. Chen and Y. Bu, ChemPhysChem, 2023, 24(11), e202200923.
16. X. Chen and Y. Bu, J. Phys. Chem. A, 2023, 127(46), 9672–9683.
17. R. E. Moore, Brit. UK Pat. Appl., Suntech. Inc., 1982, 6B2079273.
18. J. L. Adcock and M. L. Robin, J. Org. Chem., 1983, 48(18), 3128–3130.
19. X. Li, S. Li, J. Lu, H. Ren, M. Zhang and W. Zhang, Phys. Chem. Chem. Phys., 2023, 25, 15897–15904.
20. G. F. Martins, T. S. Castro and D. A. C. Ferreira, J. Mol. Model., 2023, 29, 319.
21. X. Gu and P. Jin, Inorg. Chem. Front., 2023, 10, 5634–5648.
22. M. P. Krafft and J. G. Riess, Angew. Chem., Int. Ed., 2023, 62(37), e202302942.
23. C. Titeca, T.-C. Jagau and F. De Proft, J. Chem. Phys., 2024, 160(6), 064115.
24. B. Nestmann and S. D. Peyerimhoff, J. Phys. B: At. Mol. Phys., 1985, 18(4), 615–626.
25. L. F. Williams, M. B. Yim and D. E. Wood, J. Am. Chem. Soc., 1973, 95(19), 6475–6477.
26. L. Shchegoleva, I. Bilkis and P. Schastnev, Chem. Phys., 1983, 82(3), 343–353.
27. X.-J. Hou and M.-B. Huang, J. Mol. Struct. THEOCHEM, 2003, 638(1), 209–214.
28. G. P. Paran, C. Utku and T.-C. Jagau, Phys. Chem. Chem. Phys., 2024, 26(3), 1809–1818.
29. G. Thiam and F. Rabilloud, J. Chem. Theory Comput., 2023, 19(10), 2842–2849.
30. C. Titeca, F. De Proft and T.-C. Jagau, J. Chem. Phys., 2022, 157(21), 214106.
31. N. Mardirossian and M. Head-Gordon, J. Chem. Phys., 2016, 144(21), 214110.
32. T. H. Dunning, J. Chem. Phys., 1989, 90(2), 1007–1023.
33. D. E. Woon and T. H. Dunning, J. Chem. Phys., 1993, 98(2), 1358–1371.
34. C. C. J. Roothaan, Rev. Mod. Phys., 1951, 23, 69–89.
35. G. Berthier, C. R. Hebd. Seances Acad. Sci., Ser. D, 1954, 238, 91–93.
36. J. A. Pople and R. K. Nesbet, J. Chem. Phys., 1954, 22(3), 571–572.
37. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
38. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120(1), 215–241.
39. O. Christiansen, H. Koch and P. Jørgensen, Chem. Phys. Lett., 1995, 243(5–6), 409–418.
40. Y. Jung, R. C. Lochan, A. D. Dutoi and M. Head-Gordon, J. Chem. Phys., 2004, 121(20), 9793–9802.
41. S. Grimme, J. Chem. Phys., 2003, 118(20), 9095–9102.
42. J. F. Stanton and R. J. Bartlett, J. Chem. Phys., 1993, 98(9), 7029–7039.
43. D. Rappoport, N. R. M. Crawford, F. Furche and K. Burke, Encyclopedia of Inorganic Chemistry, John Wiley & Sons, Ltd, 2009.
44. J. L. Borioni, M. Puiatti, D. M. A. Vera and A. B. Pierini, Phys. Chem. Chem. Phys., 2017, 19, 9189–9198.
45. G. Santra and J. M. L. Martin, AIP Conf. Proc., 2019, 2186(1), 030004.
46. Y. Zhang, P. M. Weber and H. Jónsson, J. Phys. Chem. Lett., 2016, 7(11), 2068–2073.
47. L. N. Anderson, M. B. Oviedo and B. M. Wong, J. Chem. Theory Comput., 2017, 13(4), 1656–1666.
48. S. Dasgupta and J. M. Herbert, J. Comput. Chem., 2017, 38(12), 869–882.
49. N. Moiseyev, Non-Hermitian Quantum Mechanics, Cambridge University Press, 2011.
50. M. Nooijen and R. J. Bartlett, J. Chem. Phys., 1995, 102(9), 3629–3647.
51. E. Epifanovsky, A. T. B. Gilbert, X. Feng, J. Lee, Y. Mao, N. Mardirossian, P. Pokhilko, A. F. White, M. P. Coons, A. L. Dempwolff, Z. Gan, D. Hait, P. R. Horn, L. D. Jacobson, I. Kaliman, J. Kussmann, A. W. Lange, K. U. Lao, D. S. Levine, J. Liu, S. C. McKenzie, A. F. Morrison, K. D. Nanda, F. Plasser, D. R. Rehn, M. L. Vidal, Z.-Q. You, Y. Zhu, B. Alam, B. J. Albrecht, A. Aldossary, E. Alguire, J. H. Andersen, V. Athavale, D. Barton, K. Begam, A. Behn, N. Bellonzi, Y. A. Bernard, E. J. Berquist, H. G. A. Burton, A. Carreras, K. Carter-Fenk, R. Chakraborty, A. D. Chien, K. D. Closser, V. Cofer-Shabica, S. Dasgupta, M. de Wergifosse, J. Deng, M. Diedenhofen, H. Do, S. Ehlert, P.-T. Fang, S. Fatehi, Q. Feng, T. Friedhoff, J. Gayvert, Q. Ge, G. Gidofalvi, M. Goldey, J. Gomes, C. E. González-Espinoza, S. Gulania, A. O. Gunina, M. W. D. Hanson-Heine, P. H. P. Harbach, A. Hauser, M. F. Herbst, M. Hernández Vera, M. Hodecker, Z. C. Holden, S. Houck, X. Huang, K. Hui, B. C. Huynh, M. Ivanov, Á. Jász, H. Ji, H. Jiang, B. Kaduk, S. Kähler, K. Khistyaev, J. Kim, G. Kis, P. Klunzinger, Z. Koczor-Benda, J. H. Koh, D. Kosenkov, L. Koulias, T. Kowalczyk, C. M. Krauter, K. Kue, A. Kunitsa, T. Kus, I. Ladjánszki, A. Landau, K. V. Lawler, D. Lefrancois, S. Lehtola, R. R. Li, Y.-P. Li, J. Liang, M. Liebenthal, H.-H. Lin, Y.-S. Lin, F. Liu, K.-Y. Liu, M. Loipersberger, A. Luenser, A. Manjanath, P. Manohar, E. Mansoor, S. F. Manzer, S.-P. Mao, A. V. Marenich, T. Markovich, S. Mason, S. A. Maurer, P. F. McLaughlin, M. F. S. J. Menger, J.-M. Mewes, S. A. Mewes, P. Morgante, J. W. Mullinax, K. J. Oosterbaan, G. Paran, A. C. Paul, S. K. Paul, F. Pavošević, Z. Pei, S. Prager, E. I. Proynov, Á. Rák, E. Ramos-Cordoba, B. Rana, A. E. Rask, A. Rettig, R. M. Richard, F. Rob, E. Rossomme, T. Scheele, M. Scheurer, M. Schneider, N. Sergueev, S. M. Sharada, W. Skomorowski, D. W. Small, C. J. Stein, Y.-C. Su, E. J. Sundstrom, Z. Tao, J. Thirman, G. J. Tornai, T. Tsuchimochi, N. M. Tubman, S. P. Veccham, O. Vydrov, J. Wenzel, J. Witte, A. Yamada, K. Yao, S. Yeganeh, S. R. Yost, A. Zech, I. Y. Zhang, X. Zhang, Y. Zhang, D. Zuev, A. Aspuru-Guzik, A. T. Bell, N. A. Besley, K. B. Bravaya, B. R. Brooks, D. Casanova, J.-D. Chai, S. Coriani, C. J. Cramer, G. Cserey, A. E. DePrince, R. A. DiStasio, A. Dreuw, B. D. Dunietz, T. R. Furlani, W. A. Goddard, S. Hammes-Schiffer, T. Head-Gordon, W. J. Hehre, C.-P. Hsu, T.-C. Jagau, Y. Jung, A. Klamt, J. Kong, D. S. Lambrecht, W. Liang, N. J. Mayhall, C. W. McCurdy, J. B. Neaton, C. Ochsenfeld, J. A. Parkhill, R. Peverati, V. A. Rassolov, Y. Shao, L. V. Slipchenko, T. Stauch, R. P. Steele, J. E. Subotnik, A. J. W. Thom, A. Tkatchenko, D. G. Truhlar, T. Van Voorhis, T. A. Wesolowski, K. B. Whaley, H. L. Woodcock, P. M. Zimmerman, S. Faraji, P. M. W. Gill, M. Head-Gordon, J. M. Herbert and A. I. Krylov, J. Chem. Phys., 2021, 155(8), 084801.
52. D. A. Matthews, L. Cheng, M. E. Harding, F. Lipparini, S. Stopkowicz, T.-C. Jagau, P. G. Szalay, J. Gauss and J. F. Stanton, J. Chem. Phys., 2020, 152(21), 214108.
53. T. A. Keith, AIMAll (Version 19.10.12), TK Gristmill Software, Overland Park KS, USA, 2019, https://aim.tkgristmill.com.
54. A. D. Becke and K. E. Edgecombe, J. Chem. Phys., 1990, 92(9), 5397–5403.
55. T. Lu and F.-W. Chen, Acta Phys.-Chim. Sin., 2011, 27(12), 2786.
56. T. Lu and F. Chen, J. Comput. Chem., 2012, 33(5), 580–592.
57. R. G. Parr and W. Yang, J. Am. Chem. Soc., 1984, 106(14), 4049–4050.
58. M. C. Oana and A. I. Krylov, J. Chem. Phys., 2007, 127(23), 124114.
59. T.-C. Jagau and A. I. Krylov, J. Chem. Phys., 2016, 144(5), 054113.
60. P. Beaujean and B. Champagne, Inorg. Chem., 2022, 61(4), 1928–1940.
61. Z. Benda, K. Rickmeyer and T.-C. Jagau, J. Chem. Theory Comput., 2018, 14(7), 3468–3478.
62. A. Shaalan Alag, D. P. Jelenfi, A. Tajti and P. G. Szalay, J. Chem. Theory Comput., 2022, 18(11), 6794–6801.
63. F. Wang, Z. Hu, X.-B. Wang, Z. Sun and H. Sun, Comput. Theor. Chem., 2021, 1202, 113295.
64. K. Sudlow and A. Woolf, J. Fluorine Chem., 1995, 71(1), 31–37.
65. K. B. Wiberg and M. Marquez, J. Am. Chem. Soc., 1998, 120(12), 2932–2938.
66. J. Jia, C. Liu, H.-S. Wu, P. V. R. Schleyer and H. Jiao, J. Phys. Chem. C, 2009, 113(19), 8077–8084.
67. P. E. Eaton and T. W. Cole, J. Am. Chem. Soc., 1964, 86(5), 962–964.
68. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14(1), 33–38.
69. R. Dennington, T. A. Keith and J. M. Millam, Gaussview Version 6, 2019.

---