Unraveling actinide–actinide bonding in fullerene cages: a DFT versus ab initio methodological study

Abstract

Actinide–actinide bonding poses a challenge for both experimental and theoretical chemists because of both the scarcity of experimental data and the exotic nature of actinide bonding due to the involvement and mixing of actinide 7s-, 6p-, 6d-, and particularly 5f-orbitals. Only a few experimental examples of An–An bonding have been reported so far. Here, we perform a methodological study of actinide–actinide bonding on experimentally known Th₂@C₈₀ and U₂@C₈₀ systems. We compared selected GGA, meta-GGA, hybrid-GGA and range-separated hybrid-GGA functionals with the results obtained using a multireference CASPT2 method, which we consider as a reference point. We show that functionals such as BP86, PBE or TPSS perform well for predicting geometries, while range-separated hybrids are superior in the description of the chemical bonding. None of the tested functionals were deemed reliable regarding the correct electronic spin ground state. Based on the results of this methodological study, we re-evaluate selected previously studied diactinide fullerene systems using more reliable protocol.

---

1. Introduction

The experimental detection of the U₂ molecule by Gorokhov et al. in 1974 started the era of actinide–actinide bonding research.¹ Only a few more systems containing An–An bonds have been reported to date. The gas-phase Th₂ was detected by means of 2D fluorescence spectroscopy.² The U(mu-H₂)UH and H₂U–UH₂ systems have been proposed to be formed in a solid argon matrix based on comparison of experimental and calculated infrared signals.³ Recent observation of a tri-thorium cluster with σ-aromatic metal–metal–metal bonding presents an example of a supported Th–Th bonding.⁴ Mass-spectroscopic evidence for various diactinide fullerenes was presented by Guo et al. in 1992 and a decade later by Akiyama et al,^5–8 but a real breakthrough was made with the full experimental characterization of endohedral fullerenes U₂@I_h(7)-C₈₀ and Th₂@I_h(7)-C₈₀ in 2018 and 2021.^9–11 Theoretical studies predicted the An–An bonds in fullerenes before they were confirmed experimentally.^12–14 Just recently, during the wrapping up of this study, Poblet and coworkers provided mass-spectroscopic evidence for a series of U₂@C_2n (2n > 50) systems prepared by laser ablation.¹⁵ They also performed theoretical calculations and revealed that U–U bonding in the studied fullerenes essentially exists in two forms, one of them being a weak U–U bond with bond order <1, r_U–U ∼ 4.0 Å, and 3 unpaired electrons at each uranium atom (septet) and the other being a triple U–U bond with r_U–U ∼ 2.5 Å and a singlet electronic state configuration, though energetically higher. Considering the abovementioned examples, actinide–actinide bonding is still rare in comparison with hundreds of existing examples of metal–metal bonds in transition metal compounds.^16,17 The fact that the experimental actinide science is more complicated is not the only reason for the scarcity of An–An bonds. Probably, the biggest obstacle is that actinides prefer binding to the ligands rather than among themselves.^18,19

The exotic nature of actinide–actinide bonding and rather poor experimental evidence of compounds containing An–An bonds inspired a number of computational studies, starting with actinide dimers^20–24 through di-actinide coordination complexes^19,25–29 to di-actinide fullerenes.^12–14,30 Effective bond orders (EBO, see Computational methods) in diatomic actinides Ac₂–Np₂ were predicted at the CASSCF and CASPT2 levels with perturbatively included spin–orbit (SO) coupling as 2, 4, and 5 in Ac₂, Th₂, and Pa₂, respectively.^21,23 Four-component calculations of U₂ with variational SO coupling predicted the EBO of four,²⁰ while the former calculations with perturbative treatment of SO coupling predicted the EBO of five.²¹ Inclusion of spin–orbit coupling was also found to slightly weaken the U–U interaction in U₂@I_h(7)-C₈₀.¹⁵

The abovementioned experimental studies of U₂@I_h(7)-C₈₀ and Th₂@I_h(7)-C₈₀ also included theoretical calculations to provide an insight into the An–An bonding in these systems. Interestingly, the CASSCF(6,6) U–U effective bond order in U₂@C₈₀ was predicted to be 0.1, while a value of 1.0 was predicted for Th–Th bonds in the Th₂@C₈₀ system.^9,10 The latter system was found to feature a classical two-electron single σ-bond (utilizing the 6s and 5d Th orbitals) with a shared electron pair in a closed-shell singlet electronic configuration. Considerably lower An–An bond order in An₂ fullerene systems as compared to the An₂ molecules is due to electron transfer from the metal to the cage. The electronegative fullerene cage acts as an electron-withdrawing ligand. For the C₈₀ cage, the experiment (XAS) shows that the oxidation state of the actinide atom enclosed is +III, while the cage bears six additional electrons.^9,10 Therefore, only some metal electrons are left for forming An–An bonds in the fullerene cage, in particular two electrons in Th₂@C₈₀ and six electrons in U₂@C₈₀.

In our previous studies of actinide fullerenes,^14,30,31 we used the DFT BP86 functional together with 60-electron core MWB potentials and a TZP-quality basis set for actinide atoms, and used the def2-SVP basis set on carbon atoms to study An–An bonding in di-actinide fullerenes. The calculations on various di-uranium fullerenes (C₆₀, C₈₀, C₈₆, C₉₀) have shown that the U–U bonding is driven by the size and shape of the fullerene cage that strongly dictates the U–U distance that in turn correlates with the U–U bond-order.¹⁴ In a follow-up, we have studied An–An bonding in An₂@C_n (An = Ac–Cm, n = 70, 80, 90) systems. In both of these studies, we used the GGA BP86 functional, as it seemed to be sufficient for the correct description of the metal–metal bonding.^13,32,33 However, there is a striking discrepancy between the CASSCF-calculated U–U EBO in U₂@C₈₀ of 0.1^9,10 and the BP86-calculated delocalization index of 1.0.^14,30 This discrepancy brought us to question how accurate are the DFT-based predictions of actinide–actinide bonds in di-actinide fullerenes.

Herein, we present a methodological study for the theoretical description of actinide–actinide bonding in actinide fullerenes. The study was performed on two experimentally characterized systems featuring An–An interactions/bonding inside a fullerene cage, U₂@I_h(7)-C₈₀ and Th₂@I_h(7)-C₈₀, vide supra. The parameters we used for the calibration method were molecular structures, interaction energies between the metal atoms and the cage, spin-state energetics, and chemical bonding indices, such as the delocalization index (DI), effective bond order (EBO), and Mayer's bond order (MBO) as well as fuzzy bond order (FBO) indicators. Of these parameters, only the molecular structures are known from the experiment. Chemical bonding is a theoretical concept and cannot be directly measured. Thus, for chemical bonding and interaction energies, we calibrated density functional theory (DFT) approaches against the best available ab initio multireference calculations, here CASPT2. We show that DFT can provide reasonable structures and energetics, fails in predicting the energy ordering of spin states, and has to be used with caution for chemical bonding parameters. Based on the results of calibration, we further revisit our former predictions of bonding along the An₂@C_n (An = Ac–Pu, n = 70, 80, 90) series of compounds³⁰ and provide new, more reliable, predictions of An–An bonding in these compounds.

2. Computational methods

Density functional theory (DFT) calculations were performed using the Gaussian 16, revision C.01 computational package.³⁴ Selection of pure DFT functionals included generalized gradient approximation (GGA) of BP86,^35,36 BLYP,^35,37 and PBE³⁸ functionals, as well as a meta-GGA TPSS³⁹ functional. Hybrid versions of functionals included B3P86,^36,40 B3LYP,^37,40 PBE0,⁴¹ and TPSSh functionals.⁴² We also evaluated the performance of range-separated CAM-B3LYP⁴³ and LC-ωHPBE functionals.⁴⁴ Local-hybrids were represented by LH14t-calPBE⁴⁵ functional as implemented in Turbomole 7.6.⁴⁶ Our attempts to include ωB97X⁴⁷ and MN15^48,49 functionals in the set failed as we were not able to achieve converged results in a reasonable time frame. DFT calculations used Grimme's empirical dispersion correction (D3) unless otherwise said.^50,51 60-electron core multireference Wood-Boring (MWB)⁵² and multireference Dirac–Fock (MDF)^53,54 pseudopotentials were used together with the corresponding basis sets that are derived from ANO-RCC basis sets for actinide atoms and correspond to (14s13p10d8f6g)/[10s9p5d4f3g]⁵⁵ and (14s13p10d8f6g)/[6s6p5d4f3g] contractions⁵⁴ for MWB60 and MDF60, respectively. The ANO-RCC-based basis sets are known to be sufficiently converged for many properties, both in the ab initio and the DFT regions, so we did not investigate the basis set convergence for An atoms.^21,56 The def2-SVP, def2-TZVP and def2-QZVP basis sets were tested for the cage carbon atoms.⁵⁷ A broken symmetry (BS) DFT approach was used to evaluate the open-shell singlet states.^58,59 In this approach, the high-spin state of the di-metal molecule of interest is calculated and then electrons on one of the metal atoms are flipped to yield an open-shell singlet, thus breaking the symmetry of the electronic state of the system.^58,59

Complete active space self-consistent field (CASSCF), complete active space 2n-order perturbation theory (CASPT2), and multi-configuration pair-density functional theory (MC-pDFT) calculations were performed using the Molcas software package, version 8.4.⁶⁰ A standard IPEA shift of 0.25 a.u. and an additional imaginary shift of 0.2 au. were used to avoid intruder states in the CASPT2 calculations. The MC-ftPBE functional was used to perform the MC-pDFT calculations. This functional has been shown to perform well for Re–Re and Cu–Cu bimetallic complexes.^61,62 The Douglas–Kroll–Hess (DKH) approximation was used together with medium Cholesky decomposition of 2-electron integrals, the ANO-RCC-VTZP (with h-type functions removed) basis set for actinide atoms, and the ANO-RCC-VDZP basis set for carbon atoms.^56,63 The basis set used for the carbon atoms presents a feasible compromise between the performance and the cost and was also used in the previously reported calculations of the EMF systems.^9,10,64 An active space of 6 electrons in 14 orbitals (6,14) was used for the U₂@C₈₀ system, while the (2,6) active space proved to be sufficient for the Th₂@C₈₀ system. Neither the larger active space nor the higher number of active electrons changed the results significantly. For more details on multireference calculations, see the ESI.†

For diactinide fullerenes, performing spin–orbit (SO) calculations presents challenges beyond the capabilities of ab initio methods. Additionally, bonding analysis using SO-DFT methods is not accessible. Previous research has suggested that SO coupling weakens the An–An bonding in U₂-based systems based on the effective bond order calculations. In Fig. S3 (ESI†), we illustrate using a model system (U₂⁶⁺) that its various spin states remain nearly degenerate, within a 3 kcal mol⁻¹ range at the scalar-relativistic Douglas–Kroll–Hess (DKH) level. When considering spin–orbit coupling, the singlet state remains the ground state, with triplet, quintet, and septet states lying at energy levels 2.6, 9.2, and 19.6 kcal mol⁻¹ above the singlet state, respectively (Fig. S3, ESI†).

Interaction energies, E_int, were calculated from the electronic energies according to eqn (1):

E_int = E(An₂@C₈₀) − (2E(An) + E(C₈₀)) (1)

The actinide atoms were considered in their ground state, i.e., [Rn]7s²6d² for thorium and [Rn]7s²5f³6d¹ for uranium.

QTAIM⁶⁵ calculations were performed using the AIMAll⁶⁶ software with both the DFT orbitals and the CASSCF/CASPT2 wavefunctions. The QTAIM parameter used to describe the strength of the interaction used in the following text is the delocalization index (DI). The DI is obtained by integrating the electron density and corresponds to a number of electron pairs shared between two atoms. Multireference wavefunctions were converted into WFX files needed for the AIMAll program using the Molden2AIM software.

Effective bond orders (EBOs) were calculated from the occupancy of the bonding and antibonding natural orbitals in the active space of the CASSCF/CASPT2 wavefunction. Like DI, EBO serves as a measure of number of electron pairs shared between the two actinides.⁶⁷ Mayer's bond orders (MBOs) and fuzzy bond orders (FBOs) were calculated using the Multiwfn software.^68–70 Natural bonding orbital (NBO) analysis was performed using the NBO 7.0 software.⁷¹ Adaptive natural density partitioning (AdNDP) analysis was performed with the AdNDP 2.0 code^72,73 using the canonical DFT as well as NBO orbitals.

3. Results and discussion

3.1. The basis set, ECP, and dispersion correction effects on geometry and interaction energy in U₂@I_h(7)-C₈₀ and Th₂@I_h(7)-C₈₀

As a first step, we investigated the effect of ECPs and basis sets on the geometry and interaction energy (E_int, eqn (1)) of the studied Th₂@C₈₀ and U₂@C₈₀ molecules. We compared older MWB ECPs with the more recent MDF ones^52–54 using the corresponding basis sets (see Computational methods). The MWB pseudopotentials are based on two-component Wood-Boring approximation, while MDF ECPs are based on the four-component Dirac–Hartree–Fock method, which should be in principle more accurate. Convergence of the basis set used for carbon atoms was studied along the def2-SVP, def2-TZVP, def2-QZVP series. Note that Th₂@C₈₀ was calculated as a singlet, although BP86 predicts a triplet ground state. This choice was done because the singlet has the lowest spin-state of Th₂@C₈₀ found by multireference calculations (Section 3.2).

The results for MWB ECPs are collected in Table 1. The calculated An–An distances in both compounds are in very good agreement with the experiment, provided the upper boundary of the U₂@C₈₀ experiment is taken (the upper boundary was recommended in the pertinent experimental work). Both the r_Th–Th and r_U–U slightly shorten with increasing size of the carbon basis set, and both are well-converged with the def2-TZVP basis set on carbon. The MDF pseudopotentials and the corresponding basis sets provide similar An–An distances and dependence on the carbon basis set as the MWB ones, as shown in Table 1. The MWB ECP gives slightly longer U–U distance than MDF ECP (3.80 vs. 3.79 Å).

Table 1 The effect of ECP on actinide and basis set on carbon on r_An−An (Å) and interaction energies E_int (kcal mol⁻¹) in Th₂@C₈₀ and U₂@C₈₀ systems. Calculated at the BP86 level

ECP	An basis	C basis	^1Th2@C80		^7U2@C80
			r	E int	r	E int
MWB60	(14s13p10d8f6g)/[10s9p5d4f3g]	def2-SVP	3.86	−377.7	3.84	−335.5
MWB60	(14s13p10d8f6g)/[10s9p5d4f3g]	def2-TZVP	3.83	−352.5	3.80	−314.3
MWB60	(14s13p10d8f6g)/[10s9p5d4f3g]	def2-QZVP	3.83	−353.9	3.80	−316.3
MDF60	(14s13p10d8f6g)/[6s6p5d4f3g]	def2-SVP	3.87	−488.4	3.82	−415.4
MDF60	(14s13p10d8f6g)/[6s6p5d4f3g]	def2-TZVP	3.84	−463.0	3.79	−394.5
MDF60	(14s13p10d8f6g)/[6s6p5d4f3g]	def2-QZVP	3.84	−464.5	3.79	−396.8
Exp.	—	—	3.82	—	3.46–3.79	—

---

Regarding the interaction energies shown in Table 1, for both kinds of ECPS and for both compounds, the def2-SVP basis set on carbon overestimates the interaction energies by about 20–25 kcal mol⁻¹, while the def2-TZVP basis set gives E_int within 2 kcal mol⁻¹ from the best def2-QZVP results. This points to converged def2-TZVP results associated with the carbon basis set and energies. In contrast, the MWB ECP-calculated E_int values are ca. 100 kcal mol⁻¹ higher than the MDF ECP values, pointing to a problem with MWB ECPs. Cross-check calculations, in which we deliberately switched the MWB and MDF basis sets revealed that MDF ECPs give consistent results with both basis sets, while MWB ECPs do not. (Table S1, ESI†). In the following discussion, we will prefer using MDF ECPs, whenever possible, specifically for Th, Pa, and U.

To see the overall effect of the used ECPs and carbon basis sets on geometry, root-mean-square deviations (RMSDs) of the calculated atomic coordinates were evaluated for both systems. The geometries are converged within RMSDs equal to or smaller than 0.002 Å for both MDF and MWB ECPs and def2-TZVP basis set (Table S2, ESI†).

The role of the dispersion interaction was investigated by performing calculations without the empirical dispersion correction (D3) for the U₂@C₈₀ system at the BP86/def2-TZVP/MDF level. Neglecting the dispersion correction lengthens the optimized U–U distance from 3.79 Å to 3.85 Å, further from the experimental value (3.79 Å) and raises the calculated interaction energy from −394.5 to −372.8 kcal mol⁻¹. The dispersion correction thus appears not to be the most crucial factor but should be definitely included in the DFT calculations of actinide fullerenes for better accuracy.

We also evaluated how ECPs and basis sets shown in Table 1 affect the selected delocalization indices (DIs). Notably, the DI is neither sensitive to the carbon basis set nor to the used ECP, as shown in Table S3 (ESI†). Based on the results in Table 1 and Tables S1–S3 (ESI†), we decided to use the more recent MDF pseudopotentials with the corresponding basis sets for actinide atoms, and the def2-TZVP basis set for carbon for further investigations. At the testing BP86 level, the chosen combination of ECP and basis sets gives well-converged An–An distances within 0.01 Å from the experiments and RMSDs vs. experimental geometries within 0.001 Å, while the interaction energies are converged within ∼2 kcal mol⁻¹ and delocalization indices within 0.01 el. pairs from the def2-QZVP results.

3.2. Electronic ground-state of Th₂@C₈₀ and U₂@C₈₀ systems. Spin-state energy differences

In systems with open-shell metal atom(s), it is important to correctly describe the correct ground state, as they are often degenerate.²³ In our previous studies, we considered the Th₂@C₈₀ system in a triplet ground state (³Th₂@C₈₀) and U₂@C₈₀ system in a septet ground state (⁷U₂@C₈₀) as predicted by the BP86 functional. Multireference calculations reveal, however, that Th₂@C₈₀ has a singlet ground state, with the lowest triplet state situated ∼10 kcal mol⁻¹ above, as shown in Table 2. This was also pointed out in a recent experimental–theoretical study, in which the lowest triplet state was found ∼ 11 kcal mol⁻¹ above the lowest singlet state at the CASPT2(2,6) level.¹⁰

Table 2 Singlet–triplet splitting in Th₂@C₈₀ and singlet–septet splitting in U₂@C₈₀ calculated at DFT, CASSCF, CASPT2 and MC-ftPBE levels using the corresponding DFT-optimized structures. Negative energy means that high-spin state (triplet for the Th₂@C₈₀ system and septet U₂@C₈₀ system) is more stable than the low-spin singlet state. Energies are in kcal mol⁻¹

	E(^3Th2@C80) − E(^1Th2@C80)				E(^7U2@C80) − E(^1U2@C80)
	DFT	CASSCF(2,6)	CASPT2	MC-ftPBE	DFT^a	CASSCF(6,14)	CASPT2	MC-ftPBE
a BS-DFT approach applied for singlet state.
BP86	−4.1	16.4	11.9	9.4	−14.4	0.0	0.3	−0.1
PBE	−2.8	17.1	12.5	10.1	−11.7	0.0	0.3	−0.1
BLYP	−3.1	15.6	11.0	8.7	−10.1	0.0	0.5	−0.2
TPSS	−4.6	16.9	12.3	9.9	−9.3	0.0	0.3	−0.1
B3P86	−3.6	15.5	11.1	8.6	−4.4	0.0	0.3	−0.1
PBE0	14.0	15.9	11.5	9.0	−2.3	0.0	0.3	−0.1
B3LYP	−1.9	15.0	10.3	8.3	−3.6	0.0	0.4	−0.1
TPSSH	−4.2	16.5	11.8	9.5	−5.8	0.0	0.3	−0.1
LC-ωHPBE	−1.7	15.7	9.9	8.4	−28.8	0.0	0.2	−0.1
CAM-B3LYP	24.7	15.1	9.6	7.9	−5.0	0.0	0.2	−0.1

---

In Table 2, we show the energy differences among the pertaining spin-states of Th₂@C₈₀ and U₂@C₈₀ molecules calculated at selected DFT and ab initio levels. The CASSCF, CASPT2, and MC-ftPBE single-point energies were calculated on top of the DFT-optimized geometries. We compare singlet–triplet splitting in Th₂@C₈₀, and singlet-septet (broken-symmetry singlet in the case of DFT) splitting in U₂@C₈₀. Note that using DFT for calculations of the triplet and quintet electronic states of U₂@C₈₀ would be unphysical, so these were not calculated at DFT levels.

For Th₂@C₈₀, all multireference methods used on any of the DFT optimized geometries predict clearly the singlet ground state that is well separated from the lowest triplet state, e.g., by 16, 12 and 9 kcal mol⁻¹ at the CASSCF, CASPT2, and MC-ftPBE levels for the BP86-optimized structure (Table 2). In contrast, most of the DFT functionals employed here predict a triplet ground state for Th₂@C₈₀. Only two of the ten functionals in Table 2 describe the spin-state ordering in Th₂@C₈₀ correctly: PBE0 (with a singlet–triplet splitting of 14.0 kcal mol⁻¹ that is comparable to the CASPT2 value of 11.5 kcal mol⁻¹) and CAM-B3LYP (with a value of 24.7 kcal mol⁻¹ as compared to the CASPT2 value of 9.6 kcal mol⁻¹). For a comparison, we also calculated the energy of the open-shell singlet for the Th₂@C₈₀ system using the broken-symmetry (BS) approach. The BS singlet lies 2.8 and 3.4 kcal mol⁻¹ above the restricted closed-shell singlet at the BP86 and PBE0 levels, respectively, and was not considered further. Note that the MC-ftPBE results are close to the CASPT2 ones, they give ca. 2 kcal mol⁻¹ smaller splitting than CASPT2.

The electronic situation is different in the U₂@C₈₀ system, where singlet and septet states are nearly degenerate, with a slight preference (<0.5 kcal mol⁻¹) for the singlet at CASSCF and CASPT2 levels and for the septet (0.1 kcal mol⁻¹) at the MC-ftPBE level. As proposed by one of the referees, this points to the importance of dynamic correlation. The CASSCF energies for the whole spin-ladder of U₂@C₈₀ are compared in Table S4 (ESI†) and are all nearly degenerate, within 1 kcal mol⁻¹. Such energy differences are on the border of accuracy to decide which spin state is the lowest one. The previously reported CASSCF(6,6) calculations predicted the singlet ground state for the U₂@C₈₀ system but details were not reported.⁹ The results in Table 2 clearly indicate that MR approaches are necessary to ensure the correct prediction of the spin state energetics in di-actinide fullerenes. Note that PBE0 is the “least wrong” functional as shown in Table 2. With certain caution, PBE0 may be used to provide qualitative estimates of the spin-state energetics in some cases, e.g., for Th₂@fullerene systems. The PBE0 functional has also proven to provide good results on energetics in fullerenes with small molecules enclosed and on spin-state energy levels in Ti@C₇₀.^74,75

Let us now have a closer look at the electronic structure of Th₂@C₈₀ and U₂@C₈₀. The CAS(2,6) wavefunction of the ¹Th₂@C₈₀ singlet is composed of three major determinants with weights of 0.89, 0.05 and 0.05, while the triplet ³Th₂@C₈₀ wavefunction is dominated by a single determinant with the weight of 0.99. The Th₂@C₈₀ system is thus to a good approximation a single-reference system. In contrast, the CAS(6,14) wavefunction of ⁷U₂@C₈₀ has one major determinant with a weight of 0.75 along with a number of minor determinants with weights below 0.01 and the singlet ¹U₂@C₈₀ wavefunction is composed of many determinants with highest weights at around 0.06. Thus, U₂@C₈₀ is a complicated multireference system. This strongly affects the bonding situation and DFT description of U–U bonding as shown in Section 3.4.

3.3. Molecular structures and interaction energies

In this Section, we further evaluate the performance of selected DFT functionals on the geometries of the ¹Th₂@C₈₀ and ⁷U₂@C₈₀ systems. DFT calculations were performed using MDF ECPs with the corresponding basis sets on an An atom and def2-TZVP basis set on cage carbons, as shown in Section 3.1. The An–An distances, the longest C–C distances, and RMSDs between structures calculated at various DFT levels and crystal structures are compared in Table 3.

Table 3 DFT-optimized and experimental r_An–An (in Å), the longest r_C–C (in Å), the difference from the experimental values (in Å, the upper value of 3.79 Å is considered for the U₂@C₈₀ system), and RMSD (in Å) between calculated and crystal structures

	^1Th2@C80					^7U2@C80
	r An–An	ΔAn–An	r C–C	ΔC–C	RMSD	r An–An	ΔAn–An	r C–C	ΔC–C	RMSD
a Range based on the results from ref. 9.
BP86	3.85	0.03	8.53	0.03	0.042	3.79	0.00	8.60	0.29	0.166
PBE	3.88	0.06	8.53	0.03	0.040	3.84	0.05	8.62	0.31	0.165
BLYP	3.81	−0.01	8.55	0.05	0.049	3.72	−0.07	8.58	0.27	0.171
TPSS	3.87	0.05	8.53	0.03	0.039	3.84	0.05	8.63	0.32	0.165
B3P86	3.87	0.05	8.46	−0.04	0.036	3.80	0.01	8.56	0.25	0.161
PBE0	3.81	−0.01	8.46	−0.04	0.036	3.85	0.06	8.57	0.26	0.162
B3LYP	3.83	0.01	8.48	−0.02	0.034	3.77	−0.02	8.56	0.25	0.164
TPSSh	3.79	−0.03	8.50	0.00	0.036	3.84	0.05	8.61	0.30	0.164
LC-ωHPBE	3.81	−0.01	8.44	−0.06	0.043	3.86	0.07	8.57	0.26	0.158
CAM-B3LYP	3.79	−0.03	8.45	−0.05	0.036	3.81	0.02	8.55	0.24	0.156
Crystal	3.82		8.50		—	3.46–3.79^a		8.31		—

---

Each of the studied systems shows qualitatively different behavior, both regarding the dependence on the DFT functional and comparison to experimental results. The DFT-optimized Th–Th distance in Th₂@C₈₀ varies between 3.79 Å (TPSSh and CAM-B3LYP) and 3.88 Å (PBE), as shown in Table 3. Notably, both the exchange and the correlation parts of a functional have an effect here. The amount of exact-exchange admixture actually seems to shorten the Th–Th distance; compare, e.g., PBE and PBE0 or TPSS with TPSSh in Table 3. But this is not true for the BLYP/B3LYP case. The effect of correlation part of the functional was observed on BP86, PBE, and BLYP levels as shown in Table 3. Range-separated hybrid functionals slightly shorten the Th–Th distance. Overall, all tested functionals show very good agreement with the experimental Th–Th distance of 3.82 Å,¹⁰ the errors are within 0.05 Å. The calculated longest C–C distance varies between 8.44 Å (LC-ωHPBE) and 8.53 Å (BP86, PBE, TPSS), again showing very good agreement with the experimental value of 8.50 Å. Overall, the Th₂@C₈₀ calculated structures differ only negligibly from the experimental data as determined from the RMSD values in Table 3.

The DFT-optimized U–U distance in the ⁷U₂@C₈₀ system varies between 3.72 (BLYP) and 3.86 Å (LC-ωHPBE). Notably, unlike Th₂@C₈₀, the exact-exchange admixture slightly elongates the U–U bond as do the range-separated hybrid LC-ωHPBE and CAM-B3LYP functionals. Overall, there are no clear trends here. The longest values of U–U bonds are predicted using pure PBE, meta-GGA TPSS, hybrid PBE0 and range-separated variant of PBE, the LC-ωHPBE. Clearly, the correlation part of the functional is more important here.

The comparison with the experiment is not simple for U₂@C₈₀. The experimental r_An–An value from XRD analysis varies between 3.46 and 3.79 Å due to several disordered X-ray structures found for the U₂@C₈₀ system, although the authors considered the shorter end of this range unreasonable.⁹ Recent theoretical calculations from the same team, however, indicated the existence of another U₂@C₈₀ minimum on a singlet hypersurface, in which the two uranium atoms form a triple bond at r_U–U ∼2.5 Å.¹⁵ This minimum lies about 10 kcal mol⁻¹ (PBE0/TZP) above the septet minimum with r_U–U ∼3.8 Å studied here. In our comparisons with the experimental data, we consider the upper experimental limit of 3.79 Å. Based on this fact, BP86, B3P86, and CAM-B3LYP give values closest to the experimental measurement. The C–C distance varies between 8.55 and 8.63 Å, which is ∼0.3 Å longer than the experimental distance of 8.31 Å. Both the U–U and C–C calculated distances are thus larger than those calculated from X-ray data for U₂@C₈₀. This is also reflected in overall larger RMSD values of ca. 0.16 Å, as shown in Table 3.

To further evaluate the performance of the DFT functionals on the structure, we calculated CASSCF, CASPT2, and MC-ftPBE single-point energies on top of the DFT-optimized geometries to determine, which DFT functional gives the lowest energy on the MR-calculated potential energy surface of Th₂@C₈₀ and U₂@C₈₀, as shown in Fig. 1. According to these results, PBE and TPSS optimized geometries give the lowest multireference (CASPT2 and MC-ftPBE) energies for both ¹Th₂@C₈₀ and ⁷U₂@C₈₀, although TPSSh, B3LYP and BP86 are also very close (up to 3 kcal mol⁻¹ for the TPSS results). Notably, here, both molecules show similar qualitative behavior.

Fig. 1 Relative single-point multireference energies (in kcal mol⁻¹) of the Th₂@C₈₀ (a) and U₂@C₈₀ (b) systems at CASSCF, CASPT2, and MC-ftPBE levels on top of the DFT-optimized geometries using CAS(2,6) and CAS(6,14) correlation spaces for Th₂@C₈₀ and U₂@C₈₀, respectively. The asterisks denote the lowest relative energy values.

The calculated interaction energies (E_int, eqn (1)) are compared in Fig. 2, and the numerical values are given in Table S5 (ESI†). Most importantly, all of the used functionals show strong metal-cage interactions in both studied systems, though the predicted values differ by tens to hundreds of kcal mol⁻¹ among the tested DFT functionals. The multireference E_int values using BP86-optimized geometries were also calculated. The CASSCF values of −274.7 and −197.7 kcal mol⁻¹ for Th₂@C₈₀ and U₂@C₈₀ are overestimated (smaller in absolute value) in comparison to the CASPT2 and MC-ftPBE values. This points to the importance of the dynamic correlation in calculations of E_int. CASPT2 results for interaction energies of −435.6 and −372.5 kcal mol⁻¹ for ¹Th₂@C₈₀ and ⁷U₂@C₈₀ are at the lower end (stronger interaction), while MC-ftPBE with −367.9 and −273.0 kcal mol⁻¹ falls between CASSCF and CASPT2. In comparison to the CASPT2 benchmark, most of the tested functionals underestimate the interaction energy but with relatively uniform trends for ¹Th₂@C₈₀, while for ⁷U₂@C₈₀, the DFT interaction energies are more spread around the CASPT2 value (Fig. 2).

Fig. 2 Comparison of the E_int (eqn (1)) values of Th₂@C₈₀ (a) and U₂@C₈₀ (b) obtained using various DFT functionals, CASSCF, CASPT2, and MC-ftPBE. The CASSCF/CASPT2/MC-ftPBE single-point calculations were performed on top of the BP86-optimized geometry with active spaces of CAS(2,6) and CAS(6,14) for Th₂@C₈₀ and U₂@C₈₀, respectively.

Except these two features, no clear trends are observed in Fig. 2. For example, for the Th₂@C₈₀ system, all but the three functionals (GGA BP86 with ca. −480 kcal mol⁻¹, hybrid B3LYP with ca. −320 kcal mol⁻¹, and RS-hybrid LC-ωHPBE with ca. −500 kcal mol⁻¹) predict comparable interaction energies (ca. −370–400 kcal mol⁻¹ while MC-ftPBE gives −370 kcal mol⁻¹) (Fig. 2a). Much larger differences among the functionals are observed for the U₂@C₈₀ system, with the lowest value of ca. −250 kcal mol⁻¹ provided by GGA BLYP and the highest value of ca. −470 kcal mol⁻¹ provided by RS-hybrid LC-ωHPBE, as shown in Fig. 2b. The performance of range-separated functionals for energetics was discussed before for energies of isomers of Sc₂@C₇₄ and Sc₂C₂@C₇₂ molecules and the importance of range-separated correction was emphasized.⁷⁶ Here, the range-separated functionals seem to overestimate E_int.

In summary, not a single functional is universal here for the structures and interaction energies. While all of the functionals have qualitatively similar performance for geometry in comparison with the experiment, the pure functionals, in particular TPSS, PBE, and BP86 seem to provide lowest energy on CASPT2 and MC-ftPBE hypersurfaces, respectively. Taking CASPT2 as a benchmark, several functionals, e.g., TPSS, PBE0, or TPSSh should provide reliable results. If a compromised functional to describe structures and interaction energies is selected here, the TPSS and PBE may be the options (Table 2, Fig. 1, and Fig. 2).

3.4. Actinide–actinide bonding

There are a number of approaches available to describe chemical bonding in molecules, although chemical bonding as such is not a directly observable quantity. Thus, we initially compared various approaches of assessing the chemical bonding, such as the QTAIM delocalization indices (DIs), Mayers bond orders (MBOs), fuzzy bond orders (FBOs), natural bonding orbitals (NBOs), and adaptive natural density partitioning (AdNDP). For details, see Computational methods. It turned out that NBO and AdNDP approaches may provide inconsistent descriptions of the An−An bonding. This is shown in ESI,† Table S7. In the following sections, we discuss DI, MBO, and FBO results.

As illustrated in Fig. 3 (numerical values in Table S8, ESI†), the DI, MBO, and FBO predict similar trends for the Th−Th bonding in Th₂@C₈₀ along the DFT axis, though the DI values are ca. 0.3 units lower than MBO and FBO values. All functionals provide qualitatively similar results here, with the trend that the Th–Th bonding slightly increases along the DFT axis as shown in Fig. 3, i.e. from pure to hybrid to range-separated functionals. Comparing the three approaches for Th₂@C₈₀, MBO and FBO values differ only marginally, while DI values are lower by ∼0.3 units. For example, at the PBE0 level, the DI gives the Th–Th bond order of ∼0.7 (i.e. the two Th atoms basically share 1.4 electrons), while MBO and FBO give values close to 1.0 and 1.1, respectively. The benchmark CASPT2(2,6) gives the DI of ∼0.71. Previous studies predicted the effective bond order of 1.0 for the Th₂@C₈₀ system using the CASSCF(2,6) wavefunction.¹⁰ Our comparative EBO calculations predict the same value (see Table S9, ESI†).

Fig. 3 Comparison of DI, MBO, and FBO values representing the Th−Th bonding in Th₂@C₈₀ calculated using the DFT, CASSCF(2,6) and CASPT2(2,6) orbitals.

The description of U–U bonding in U₂@C₈₀ (Fig. 4, numerical values in Table S8, ESI†) is qualitatively different from that of Th₂@C₈₀ (compare with Fig. 3). The DFT-predicted FBO values are systematically higher by ∼0.3 units and MBO values are systematically lower by ∼0.3 units when compared to DI. Furthermore, the trends along the DFT axis are rather opposite, and the U–U interaction decreases when going from pure through hybrid to range-separated DFT functionals.

Fig. 4 Comparison of DI, MBO, and FBO values representing the U−U interaction in U₂@C₈₀ calculated using the DFT, CASSCF(6,14) and CASPT2(6,14) orbitals.

Of all the bonding analyses, we trust most the QTAIM-based delocalization index (DI) that is based on the integration of the electron density and also accessible for MR methods. Importantly, DIs could be obtained also at ab initio multireference CASSCF and CASPT2 levels (see Computational methods). Let us now have a look at DIs in Th₂@C₈₀ (Fig. 3, numerical values in Table S8, ESI†) in detail. As indicated above, the delocalization index for the Th–Th bonding (DI_Th–Th) is relatively stable for all of the tested DFT functionals, ranging from ∼0.60 for PBE and GGA functionals to ∼0.75 for CAM-B3LYP and LC-ωHPBE functionals (Fig. 3). This is in agreement with the benchmark DI_Th–Th values, calculated at the ab initio CASSCF and CASPT2 levels of 0.73 and 0.71, respectively.

A closer look at Fig. 3 reveals that the exact-exchange admixture slightly increases the DI_Th–Th in the case of PBE-based functionals; compare PBE (0.60) vs. PBE0 (0.69), but there is no effect when comparing B and B3 including functionals; compare BLYP (0.69) vs. B3LYP (0.69). Six out of the ten tested functionals, the BLYP, PBE0, B3LYP, TPSSh, LC-ωHPBE, and CAM-B3LYP functionals with DI_Th–Th values of 0.69, 0.69, 0.69, 0.73, 0.76 and 0.75, respectively, give values that are reasonably close to the CASPT2 benchmark of 0.71. Even though the various functionals describe the Th–Th interaction in Th₂@C₈₀ with small differences, all of them consistently evaluate that there is a bonding interaction of an order of at least ∼0.6 – that means 0.6 of an electron pair is shared between the two Th atoms. For example, LC-ωHPBE gives a DI of 0.81. These results are in line with the recent experimental–theoretical study results, where the CASSCF(2,6) results concluded a single sigma bond between the two Th atoms (based on the MO picture and EBO of 0.9).¹⁰

Notably, the lowest triplet ³Th₂@C₈₀ (see above) also features Th–Th bonding interaction as predicted previously by us.³⁰ In the triplet state, the bonding is realized via two one-electron-two-center (OETC) bonds. The DFT-predicted DI_Th–Th in ³Th₂@C₈₀ is around ∼0.8 for a few tested functionals, such as BP86, TPSS, and TPSSh. For details on bonding in ³Th₂@C₈₀, see ref. 30.

The DI_Th-C80 values (Table S8, ESI†) that represent overall interaction of a single Th with the cage fall in a range of 3.2–3.75 at the DFT level; the lowest reported values (∼3.2) are obtained using the range-separated and hybrid DFT, while GGA functionals give the highest values (∼3.7). Approximately 50% of the Th-C₈₀ DI value comes from the interaction between the thorium atom and the six closest carbon atoms. The DI_Th-C80 values obtained from CASSCF(2,6) and CASPT2(2,6) wavefunctions are ca. one unit lower, 2.81 and 2.65, respectively. This may in part originate from a smaller DZ-quality basis set on carbon but using a larger basis set fullerene cage would be computationally prohibitive in MR calculations. Note that while larger DI_Th–Th values are obtained using range-separated functionals, the opposite happens for the DI_Th-C80, i.e., the stronger the An–An interaction, the weaker the interaction with the cage (Fig. 4 and Table S8, ESI†).

There is a mild correlation between the DI_Th–Th and the Th–Th distance calculated with different functionals, with R² of 0.83 (Fig. 5a). The shorter the Th–Th distance, the stronger the bonding interaction (expressed as DI_Th–Th), even though different DFT functionals are used. Similar correlation has been previously studied with the BP86 functional in connection with the DI_U–U interaction and cage size. According to our previous study, a larger cage implies a longer An–An distance, which in turn results in a weaker An−An interaction (lower DI_An–An).¹⁴ An opposite trend – with the R² of 0.74 – is seen when correlating the DI_Th-C80 and the distance between two farthest carbon atoms (Fig. 5b).

Fig. 5 Correlation between the DI_Th–Th and the Th–Th distance (a) and between the DI_Th-C80 and the longest distance between two carbon atoms (b), based on the results calculated using different functionals.

The description of bonding in the U₂@C₈₀ system is more complicated in comparison to that in Th₂@C₈₀. While the span of r_U–U values in U₂@C₈₀ is still relatively small, ∼0.14 Å (Table 2), the calculated DI_U–U strongly varies with the DFT functionals, starting from ∼1.1 for GGA functionals all the way down to 0.17 for the range-separated LC-ωHPBE level (Fig. 4). The benchmark DI_U–U values calculated at the multireference CASSCF and CASPT2 levels of 0.07 for both show only a marginal interaction. The DI_U–U of ∼1.1 for GGA functionals and ∼0.87 for meta-GGA TPSS propose an equivalent of a single bond as observed in the previous work, where the DI_U–U of 1.0 was predicted at the BP86/MWB/def2SVP level and geometry.^14,30 However, the exact-exchange admixture in hybrid functionals brings the DI_U–U down to 0.62 for TPSSh (10% HF exchange) and 0.36 for PBE0 (25% HF exchange). The range-separated functionals further lower the DI to 0.23 at the CAM-B3LYP level and 0.17 for the LC-ωHPBE functional. Thus, only the range-separated functionals seem to provide the bonding picture in U₂@C₈₀ close to the ab initio benchmark, though the interaction is still overestimated (DI_U–U of 0.07 vs. ∼0.2).

Regarding DI_U–C80, approximately 50% of this value arises from the interaction between uranium atoms and the six closest carbon atoms, as observed in Th₂@C₈₀. The DI_U–C80 values obtained from CASSCF and CASPT2 wavefunctions are 2.64 and 2.46, respectively (Table S8, ESI†). This is 1 to 1.5 lower than the DFT predictions. All functionals overestimate this value, with the GGA functionals doing so to a greater extent. For example, BP86 gives a DI of 4.05 and CAM-B3LYP gives a DI of 3.26. Nevertheless, the progression of DI_U-C80 towards the range-separated functionals shows that not only U–U interaction but also U–C interaction is overestimated by the pure GGA functionals and this decreases with the inclusion of exact exchange in hybrid functionals and even more so for range-separated hybrids.

The correlation between the calculated r_Th–Th and DI_Th–Th observed for Th₂C₈₀ (Fig. 5a) is not repeated for U₂@C₈₀ (Fig. 6a). The strength of the interaction as described using the DI appears to be dependent on the functional used, rather than on r_U–U. This further supports the conclusion that DFT has a problem with the description of the U–U bonding in U₂@C₈₀. We speculate that this system converges to different microstates but this is difficult to prove. Mild correlation (R² = 0.65) is seen between the DI_U-C80 and the longest C–C distance, similar to the case of Th₂@C₈₀. The role of exact exchange was further evaluated. The BP86 functional was modified with growing percentage of exact exchange increments starting at 10% and ending at 90%. The higher the portion of exact exchange in the functional, the higher the bond orders in the Th₂C₈₀ molecule. The results for U₂@C₈₀ were inconclusive (see Fig. S4 and commentary in the ESI†).

Fig. 6 Dependence of DI_U–U on r_U–U (a) and of the DI_U-C80 on the longest distance between two carbon atoms (b).

The dependence of QTAIM parameters on a functional in actinide compounds was discussed before. Slight U–Cl bond DI overestimation was shown while using pure functionals to calculate the Cs₂UO₂Cl₄ system.⁷⁷ During finalization of this work, we also tested a local-hybrid LHT14-calPBE functional available in Turbomole on both ¹Th₂@C₈₀ and ⁷U₂@C₈₀ systems. Subsequent QTAIM analysis gave DI_An–An values of 0.68 and 0.98 for Th₂@C₈₀ and U₂@C₈₀ systems. These results are comparable with the performance of pure functionals, which is somehow unexpected. The tested local-hybrid functional thus did not improve the description of U–U bonding.

3.5. Re-evaluation of bonding in previously reported An₂@C_n systems

With the methodology assessed, we revisited our former BP86/def2-SVP predictions and recalculated previously studied An₂@C_n systems, where possible.³⁰ The following protocols were used: first, the geometries were optimized using the TPSS-D3 functional, the def2-TZVP basis set on the carbon atoms, and MDF (Ac–U) or MWB (MDF is not available beyond U) pseudopotential with the corresponding basis set on actinide atoms. Second, the LC-ωHPBE single-point calculations were performed on TPSS-D3 geometries and the obtained electron density was analyzed using QTAIM to calculate the required DI parameters. Because even the best LC-ωHPBE functional overestimates the An–An bonding interactions (Fig. 4), we also performed QTAIM bonding analysis at the CASSCF level for comparison. Active spaces of CAS(2,6), CAS(4,14), CAS(6,14), and CAS(10,14) were used for Th₂@C_n, Pa₂@C_n, U₂@C_n, and Pu₂@C_n systems, respectively. ANO-RCC-9s8p6d4f2g basis sets were used for the actinide atoms and ANO-RCC-VDZP basis set for carbon atoms. The CASPT2 level was used to determine the electronic ground states, as shown in Table S10 (ESI†).

The results are collected in Table 4. For the closed-shell ¹Th₂@C_n systems, excellent agreement is observed between previously reported results and newly obtained DIs, both at the LC-ωHPBE and CASSCF level. The Pa–Pa interaction in ⁵Pa₂@C₈₀ that was previously reported be the strongest interaction between two actinide atoms enclosed in C₈₀, is actually comparable to the Th–Th interaction in the Th₂@C₈₀ system, with the DI_Pa-Pa of 0.8 calculated from the LC-ωHPBE orbitals and that of 0.7 calculated from the CASSCF(4,14) wavefunction.

Table 4 Comparison of previous^a and present DI_An–An values in M₂@C_n fullerenes. The electronic ground states of compounds were derived from the CASPT2 calculations

QTAIM	BP86	LC-ωHPBE-D3	CASSCF
a (BP86/def2-SVP/MWB/MDF). b Note that former BP86 calculations considered triplet ground state ^3Th2@Cn.
Geometry	BP86^a^,14,30	TPSS-D3	TPSS-D3
^1Th2@C70	0.4^b	0.5	0.6
^1Th2@C80	0.8^b	0.8	0.7
^1Th2@C90	0.0^b	0.0	0.0
^5Pa2@C80	1.3	0.8	0.7
^3U2@C60	2.1	1.6	1.2
^7U2@C70	0.7	0.4	0.1
^7U2@C80	1.0	0.3	0.1
^1U2@C90	0.1	0.0	0.0
^11Pu2@C90	0.4	0.0	0.0

---

In the U₂@C_70,80,90 systems (Table 4), septet and singlet electronic ground states are nearly degenerate. The U₂@C₇₀ and U₂@C₈₀ systems slightly prefer a septet state, while U₂@C₉₀ prefers a singlet ground state. While the DIs of 0.3–0.4 are calculated at the DFT level for the U₂@C₇₀ and U₂@C₈₀ system, analysis of the CASSCF wavefunction shows almost negligible U–U interaction of ∼0.1. Negligible U–U interaction is calculated consistently using both DFT and CASSCF for U₂@C₉₀. The Pu₂@C₉₀ system was previously thought to feature a long Pu–Pu interaction at r_Pu–Pu = 5.9 Å. Both the LC-ωHPBE and the CASSCF results show no interaction between the two plutonium atoms.

We also added the U₂@C₆₀ system in Table 4 and Fig. S5 (ESI†) as it has received considerable attention in the past. Its mass spectroscopic signal was first reported by Guo in 1992.⁵ In 2006, Wu and Lu predicted six one-electron two-center bonds between two uranium atoms in the system, based on DFT calculations.¹² In 2008, Infante et al. argued that the interaction is a mere artifact due to the small size of the cage.¹³ Here, we conclude that DFT calculations tend to overestimate the strength of U–U bonding in ³U₂@C₆₀, with the CASSCF(6,14) level predicting a DI_U–U of 1.2.

Conclusions

In this work, we studied actinide–actinide bonding in experimentally known Th₂@C₈₀ and U₂@C₈₀ systems using DFT and multireference CASSCF, CASPT2 and MC-ftPBE methods to calibrate the theoretical methodology for actinide–actinide bonding in fullerenes. Clear differences between the two studied systems were observed. The benchmark level CASPT2(2,6) predicts a closed-shell singlet electronic ground state for Th₂@C₈₀ with a two-electron single Th–Th bond and a single dominant determinant in the CASSCF wavefunction. For this system, the DFT works generally well and provides consistent results for the geometry, energetics, and various chemical bonding indicators. However, even in this seemingly single-reference molecule, most of the tested DFT functionals fail to predict the correct spin ground state in Th₂@C₈₀. All used functionals predict a triplet ground state, except for the PBE0 and CAM-B3LYP functionals that predict a closed shell singlet ground state. However, we consider this incidental.

For U₂@C₈₀, the benchmark CASPT2(6,14) calculations predict a rather degenerate system in which the three 5f electrons on each uranium can almost freely combine their spins in nearly degenerate septet, quintet, triplet, and singlet electronic states, with a slight preference for the ground state septet. The corresponding wavefunctions are composed of a number of low-weight determinants. Here, the DFT results can still be reasonable for the structure and energetics, but chemical bonding indicators, such as the delocalization index are strongly functional-dependent. Also, none of the tested functionals provides the correct spin ground state.

Overall, the pure GGA and meta-GGA functionals, such as BP86, PBE and TPSS show a better performance for the prediction of molecular structures when compared to experiment and MR calculations, while the range-separated functionals describe better the electron density – particularly the chemical bonding – more accurately, although they still overestimate the DI obtained from MR approaches. Thus, the DFT can be used here for qualitative estimates of bonding interactions. However, for obtaining a correct electronic ground state, one should, in any case, employ multireference calculations with appropriate CAS space. The DFT proved here to be unreliable, regardless of whether GGA, hybrid-GGA or range-separated hybrid-GGA functionals are used, even for the nearly single-reference Th₂@C₈₀ system.

We have also shown that the delocalization index, Mayer bond order, and Fuzzy bond order indices provide qualitatively same trends as the QTAIM delocalization index, though the absolute values are different. Thus, a low-cost analysis, such as MBO, can be used to pre-evaluate chemical bonding instead of much more expensive QTAIM.

Finally, we revisited calculations of the previously reported An₂@C_n systems using the best affordable levels acquired from the presented calibration, both DFT and ab initio (CASSCF for DI, CASPT2 for spin state ordering) to consolidate the proposed approach. We conclude that while the previously reported Th₂@C₇₀ and Pa₂@C₈₀ systems feature the actinide–actinide bonds as predicted, the U–U and Pu–Pu interactions are actually almost negligible in Pu₂@C₉₀, U₂@C₇₀, U₂@C₈₀, and U₂@C₉₀ systems. Contrary to the previous DFT predictions of six one-electron-two-center bonds U–U bonds in U₂@C₆₀, this system features an equivalent of a single U–U bond at the CASPT2 level.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank Dr. Stanislav Komorovský for fruitful discussions about spin–orbit effects. This project was supported by Czech Science Foundation grant 21-17806. Computational resources were supplied by the project “e-Infrastruktura CZ” (e-INFRA CZ LM2018140) supported by the Ministry of Education, Youth and Sports of the Czech Republic.

References

1. L. Gorokhov, A. Emelyanov and Y. S. Khodeev, High Temp., 1974, 12, 1156–1158.
2. T. Steimle, D. L. Kokkin, S. Muscarella and T. Ma, J. Phys. Chem. A, 2015, 119, 9281–9285.
3. P. F. Souter, G. P. Kushto, L. Andrews and M. Neurock, J. Am. Chem. Soc., 1997, 119, 1682–1687.
4. J. T. Boronski, J. A. Seed, D. Hunger, A. W. Woodward, J. van Slageren, A. J. Wooles, L. S. Natrajan, N. Kaltsoyannis and S. T. Liddle, Nature, 2021, 598, 72–75.
5. T. Guo, M. D. Diener, Y. Chai, M. J. Alford, R. E. Haufler, S. M. McClure, T. Ohno, J. H. Weaver, G. E. Scuseria and R. E. Smalley, Science, 1992, 257, 1661–1664.
6. K. Akiyama, Y. Zhao, K. Sueki, K. Tsukada, H. Haba, Y. Nagame, T. Kodama, S. Suzuki, T. Ohtsuki, M. Sakaguchi, K. Kikuchi, M. Katada and H. Nakahara, J. Am. Chem. Soc., 2001, 123, 181–182.
7. K. Akiyama, K. Sueki, K. Tsukada, T. Yaita, Y. Miyake, H. Haba, M. Asai, T. Kodama, K. Kikuchi, T. Ohtsuki, Y. Nagame, M. Katada and H. Nakaharaa, J. Nucl. Radiochem. Sci., 2002, 3, 151–154.
8. K. Akiyama, K. Sueki, H. Haba, K. Tsukada, M. Asai, T. Yaita, Y. Nagame, K. Kikuchi, M. Katada and H. Nakahara, J. Radioanal. Nucl. Chem., 2003, 255, 155–158.
9. X. Zhang, Y. Wang, R. Morales-Martinez, J. Zhong, C. de Graaf, A. Rodriguez-Fortea, J. M. Poblet, L. Echegoyen, L. Feng and N. Chen, J. Am. Chem. Soc., 2018, 140, 3907–3915.
10. J. Zhuang, R. Morales-Martínez, J. Zhang, Y. Wang, Y.-R. Yao, C. Pei, A. Rodríguez-Fortea, S. Wang, L. Echegoyen, C. de Graaf, J. M. Poblet and N. Chen, Nat. Commun., 2021, 12, 2372.
11. W. Cai, C.-H. Chen, N. Chen and L. Echegoyen, Acc. Chem. Res., 2019, 52, 1824–1833.
12. X. Wu and X. Lu, J. Am. Chem. Soc., 2007, 129, 2171–2177.
13. I. Infante, L. Gagliardi and G. E. Scuseria, J. Am. Chem. Soc., 2008, 130, 7459–7465.
14. C. Foroutan-Nejad, J. Vícha, R. Marek, M. Patzschke and M. Straka, Phys. Chem. Chem. Phys., 2015, 17, 24182–24192.
15. A. Moreno-Vicente, Y. Roselló, N. Chen, L. Echegoyen, P. W. Dunk, A. Rodríguez-Fortea, C. de Graaf and J. M. Poblet, J. Am. Chem. Soc., 2023, 145, 6710–6718.
16. F. A. Cotton, N. F. Curtis, C. B. Harris, B. F. G. Johnson, S. J. Lippard, J. T. Mague, W. R. Robinson and J. S. Wood, Science, 1964, 145, 1305–1307.
17. J. F. Berry and C. C. Lu, Inorg. Chem., 2017, 56, 7577–7581.
18. S. T. Liddle and D. P. Mills, Dalton Trans., 2009, 5592–5605.
19. G. Cavigliasso and N. Kaltsoyannis, Inorg. Chem., 2006, 45, 6828–6839.
20. S. Knecht, H. J. A. Jensen and T. Saue, Nat. Chem., 2019, 11, 40–44.
21. L. Gagliardi and B. O. Roos, Nature, 2005, 433, 848–851.
22. H. S. Hu and N. Kaltsoyannis, Phys. Chem. Chem. Phys., 2017, 19, 5070–5076.
23. B. O. Roos, P.-Å. Malmqvist and L. Gagliardi, J. Am. Chem. Soc., 2006, 128, 17000–17006.
24. L. Gagliardi, P. Pyykkö and B. O. Roos, Phys. Chem. Chem. Phys., 2005, 7, 2415–2417.
25. M. Straka and P. Pyykkö, J. Am. Chem. Soc., 2005, 127, 13090–13091.
26. B. O. Roos and L. Gagliardi, Inorg. Chem., 2006, 45, 803–807.
27. G. La Macchia, M. Brynda and L. Gagliardi, Angew. Chem., Int. Ed., 2006, 45, 6210–6213.
28. G. Cavigliasso and N. Kaltsoyannis, Dalton Trans., 2006, 5476–5483.
29. C. Z. Wang, J. K. Gibson, J. H. Lan, Q. Y. Wu, Y. L. Zhao, J. Li, Z. F. Chai and W. Q. Shi, Dalton Trans., 2015, 44, 17045–17053.
30. A. Jaroš, C. Foroutan-Nejad and M. Straka, Inorg. Chem., 2020, 59, 12608–12615.
31. J. Kaminský, J. Vícha, P. Bouř and M. Straka, J. Phys. Chem. A, 2017, 121, 3128–3135.
32. C. J. Barden, J. C. Rienstra-Kiracofe and H. F. S. III, J. Chem. Phys., 2000, 113, 690–700.
33. S. Yanagisawa, T. Tsuneda and K. Hirao, J. Chem. Phys., 2000, 112, 545–553.
34. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, 2016.
35. A. D. Becke, Phys. Rev. A, 1988, 38, 3098–3100.
36. J. P. Perdew, Phys. Rev. B, 1986, 33, 8822–8824.
37. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B, 1988, 37, 785–789.
38. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1997, 78, 1396.
39. J. Tao, J. P. Perdew, V. N. Staroverov and G. E. Scuseria, Phys. Rev. Lett., 2003, 91, 146401.
40. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
41. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
42. V. N. Staroverov, G. E. Scuseria, J. Tao and J. P. Perdew, J. Chem. Phys., 2003, 119, 12129–12137.
43. T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51–57.
44. T. M. Henderson, A. F. Izmaylov, G. Scalmani and G. E. Scuseria, J. Chem. Phys., 2009, 131, 044108.
45. A. Arbouznikov and M. Kaupp, J. Chem. Phys., 2014, 141, 204101.
46. TURBOMOLE V7.6, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; available from https://www.turbomole.com.
47. J.-D. Chai and M. Head-Gordon, J. Chem. Phys., 2008, 128, 084106.
48. H. S. Yu, X. He, S. L. Li and D. G. Truhlar, Chem. Sci., 2016, 7, 5032–5051.
49. H. S. Yu, X. He and D. G. Truhlar, J. Chem. Theory Comput., 2016, 12, 1280–1293.
50. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
51. L. Goerigk, A. Hansen, C. Bauer, S. Ehrlich, A. Najibi and S. Grimme, Phys. Chem. Chem. Phys., 2017, 19, 32184–32215.
52. W. Küchle, M. Dolg, H. Stoll and H. Preuss, J. Chem. Phys., 1994, 100, 7535–7542.
53. M. Dolg and X. Cao, J. Phys. Chem. A, 2009, 113, 12573–12581.
54. A. Weigand, X. Cao, T. Hangele and M. Dolg, J. Phys. Chem. A, 2014, 118, 2519–2530.
55. X. Cao, M. Dolg and H. Stoll, J. Chem. Phys., 2003, 118, 487–496.
56. B. O. Roos, R. Lindh, P.-Å. Malmqvist, V. Veryazov and P.-O. Widmark, Chem. Phys. Lett., 2005, 409, 295–299.
57. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
58. L. Noodleman, J. Chem. Phys., 1981, 74, 5737–5743.
59. F. Neese, J. Phys. Chem. Solids, 2004, 65, 781–785.
60. F. Aquilante, J. Autschbach, R. K. Carlson, L. F. Chibotaru, M. G. Delcey, L. De Vico, I. F. Galván, N. Ferré, L. M. Frutos, L. Gagliardi, M. Garavelli, A. Giussani, C. E. Hoyer, G. Li Manni, H. Lischka, D. Ma, P. Å. Malmqvist, T. Müller, A. Nenov, M. Olivucci, T. B. Pedersen, D. Peng, F. Plasser, B. Pritchard, M. Reiher, I. Rivalta, I. Schapiro, J. Segarra-Martí, M. Stenrup, D. G. Truhlar, L. Ungur, A. Valentini, S. Vancoillie, V. Veryazov, V. P. Vysotskiy, O. Weingart, F. Zapata and R. Lindh, J. Comput. Chem., 2016, 37, 506–541.
61. R. K. Carlson, D. G. Truhlar and L. Gagliardi, J. Chem. Theory Comput., 2015, 11, 4077–4085.
62. A. Stańczak, J. Chalupský, L. Rulíšek and M. Straka, Chem. Phys. Chem., 2022, 23, e202200076.
63. B. O. Roos, R. Lindh, P.-Å. Malmqvist, V. Veryazov and P.-O. Widmark, J. Phys. Chem. A, 2004, 108, 2851–2858.
64. X. Li, Y. Roselló, Y.-R. Yao, J. Zhuang, X. Zhang, A. Rodríguez-Fortea, C. de Graaf, L. Echegoyen, J. M. Poblet and N. Chen, Chem. Sci., 2021, 12, 282–292.
65. R. F. W. Bader, Chem. Rev., 1991, 91, 893–928.
66. Todd. A. Keith, AIMAll (Version 19.10.12) (version Version 19.10.12) TK Gristmill Software, Overland Park KS, USA 2019.
67. B. O. Roos, A. C. Borin and L. Gagliardi, Angew. Chem., Int. Ed., 2007, 46, 1469–1472.
68. I. Mayer, Chem. Phys. Lett., 1983, 97, 270–274.
69. I. Mayer and P. Salvador, Chem. Phys. Lett., 2004, 383, 368–375.
70. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
71. E. D. Glendening, J. K. Badenhop, A. E. Reed, J. E. Carpenter, J. A. Bohmann, C. M. Morales, P. Karafiloglou, C. R. Landis and F. Weinhold, NBO 7.0 Theoretical Chemistry Institute, University of Wisconsin, Madison 2018.
72. D. Y. Zubarev and A. I. Boldyrev, Phys. Chem. Chem. Phys., 2008, 10, 5207–5217.
73. N. V. Tkachenko and A. I. Boldyrev, Phys. Chem. Chem. Phys., 2019, 21, 9590–9596.
74. L. Tučková, A. Jaroš, C. Foroutan-Nejad and M. Straka, Phys. Chem. Chem. Phys., 2023, 25, 14245–14256.
75. A. Jaroš, M. Sasar, L. Tučková, E. F. Bonab, Z. Badri, M. Straka and C. Foroutan-Nejad, Adv. Electron. Mater., 2023, 9, 2300360.
76. R. Zhao, Y. Guo, P. Zhao, M. Ehara, S. Nagase and X. Zhao, J. Phys. Chem. C, 2016, 120, 1275–1283.
77. J. P. W. Wellington, A. Kerridge and N. Kaltsoyannis, Polyhedron, 2016, 116, 57–63.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3cp03606e

---