Single-ion magnet behaviour in highly axial lanthanide mononitrides encapsulated in boron nitride nanotubes: a quantum chemical investigation

Abstract

Lanthanide-based single-ion magnets (Ln-SIMs) have garnered significant interest for their potential application in molecular-level information storage devices. Among various strategies to enhance the magnetization blocking barrier in SIMs, synthesizing highly axially symmetric compounds is the most promising approach. In the present work, using state-of-the-art computational tools, we have thoroughly examined the electronic structure, bonding, and magnetic anisotropy of lanthanide mononitrides [LnN] (where Ln = Dy(III) and Tb(III)) and their encapsulation in zigzag boron nitride nanotubes (BNNTs) with diameters of (8,0) and (9,0) to explore novel hybrid assemblies. Using periodic density functional theory calculations and energy decomposition analysis, we have thoroughly analyzed the structural and energetic perspectives towards encapsulation of [LnN] molecules in parallel and perpendicular modes in BNNT(8,0) (8Ln_|| and 8Ln_⊥) and BNNT(9,0) (9Ln_|| and 9Ln_⊥) tubes. Binding energy calculations suggest that the parallel arrangement of [LnN] is energetically more favourable (>30 kJ mol⁻¹) than the perpendicular arrangement, with the BNNT(8,0) tube being energetically more preferred over the BNNT(9,0) tube for encapsulation. Non-covalent interaction plots clearly show dominant van der Waals interactions in 8Dy_||/8Tb_||, stabilizing it compared to other assemblies. CASSCF calculations suggest that both [TbN] and [DyN] show a pure Ising-type ground state with a giant barrier height of >1800 cm⁻¹ and strictly no ground-state quantum tunnelling of magnetization. CASSCF calculations predict that the 8Dy_|| and 8Tb_|| assemblies show record high ab initio blockade barrier (U_cal) values of ∼1707 and 1015 cm⁻¹, respectively. Although 9Dy_⊥ is an energetically unfavourable mode, this orientation benefits from the tube's crystal field, which leads to a U_cal value of ∼1939 cm⁻¹, suggesting that encapsulation could further enhance the U_cal values. Contrarily, the [TbN] molecules show a dramatic increase in the tunnel splitting values upon encapsulation in BNNT tubes, leading to a drastic decrease in U_cal values. Our in silico strategy offers insights into the magnetic anisotropy of simple [DyN] and [TbN] molecules and possible ways to integrate these molecules into BNNTs to generate hybrid magnetic materials for information storage applications.

---

Introduction

Strongly anisotropic mononuclear lanthanide complexes are regarded as the best molecule-based magnets that exhibit magnetic relaxation close to the boiling point of liquid nitrogen, opening up great potential for high-performance information storage, quantum information processing, and spintronic device applications.^1–3 Exceptionally stable trivalent states and simple structures associated with mononuclear complexes make them ideal for fine-tuning magnetic anisotropy through robust magneto-structural correlation.^4–8 Designing highly anisotropic Ln-SIMs requires a careful ligand field design depending on the electron cloud shape, where the prolate shape density requires a complete equatorial ligand field and the oblate shape density necessitates a complete axial ligand field to stabilize the largest m_J as the ground state.^4–6,9 Over the past decade, magnetochemists have significantly increased the barrier height (U_eff) and blocking temperature (T_B) of Ln(III) based SIMs, which determines the performance of an SMM.^10–13 The best-performing single-molecule magnets (SMMs) are the family of organometallic dysprosocenium complexes, where the [Dy(Cp*)(Cp^iPr5)][B(C₆F₅)₄], [Cp* denotes pentamethyl-cyclopentadienyl, Cp^ip5 = penta-iso-propylcyclopentadienyl] complex displays a giant barrier height over 1500 cm⁻¹ with a T_B value of just 80 K.¹⁰ Increasing the blocking temperature and achieving stability under ambient conditions are imperative for utilizing these molecules for any possible applications in information storage devices. In general, a large coordination number (CN > 7) is required to stabilize air-stable Dy(III) complexes¹⁴ and to achieve SIM behaviour in these complexes, the symmetry around the Dy(III) ion and the strength of the axial and equatorial ligand field must be balanced. The presence of axial symmetry in square antiprismatic/axially compressed octahedral (D_4d/D_4h),¹⁵ pentagonal bipyramidal (D_5h), and hexagonal bipyramidal (D_6h) environments with short Dy–L_ax bonds compared to the Dy–L_eq bonds is one feasible way to generate an axial ligand field in high-coordination-number complexes.^16–21 Significant research highlights include the discovery of [Dy(O^tBu)₂(L)₄]⁺, with an axially compressed octahedral (D_4h) configuration, [Dy(O^tBu)₂(py)₅]⁺, featuring a pentagonal bipyramidal (PBP) geometry (D_5h), and R/S-[Dy(L^N6)(OSiPh₃)₂], in a D_6h environment, with U_eff values over 1800 K in all the cases.^{16,19,22,23–28} Xu et al. have reported a classic example of air-stable bis-alkoxide Dy(III) crown ether complexes where a weaker ligand field offered by neutral 18-crown-6 ligands helps in generating a giant U_eff value of 2427 K.²⁹ Recently, we predicted a series of Dy(III) complexes comprised of unsaturated 18-crown-6 ligands displaying a highly symmetric D_6h geometry around the Dy(III) ion with an ab initio blockade barrier over ∼3000 K.⁸ These complexes show a promising avenue for the design of highly anisotropic, air-stable SIMs with a giant magnetic relaxation barrier height.

Besides designing highly anisotropic SMMs, there have been several efforts to understand the encapsulation of SMMs into hollow/porous materials to create hybrid systems retaining SMM properties, which can finally be integrated into devices for reading and writing purposes. Numerous Ln-based SIMs/SMMs have been studied in detail by encapsulating them in fullerenes/nanotubes or depositing them on 2D nanostructures to obtain hybrid assemblies.^30–33 In this context, 1D nanotubes are emerging materials that allow the deposition/encapsulation of SMMs on the surface of the tube/inside the tube. Tangoulis et al. reported the functionalization of carbon nanotubes (CNTs) with the [Dy₂(NO₃)₂(saph)₂(DMF)₄] complex (where H₂saph = N-salicylidene-o-aminophenol) demonstrating that the grafted molecules maintain their SMM characteristics.³⁴ In one of the works, Nakanishi et al. encapsulated DySc₂N@C₈₀-SMMs in single-walled carbon nanotubes (SWCNTs), which exhibits suppressed quantum tunnelling of magnetization (QTM).³⁵ Popov and co-workers have thoroughly explored endohedral metallofullerenes where Ln(III) atoms are clubbed with non-metals,^36,37 resulting in giant magnetic anisotropy and magnetic exchange coupling, an essential criteria for building high-performance SMMs. Using a capillary method, Yamashita and co-workers have recently demonstrated the encapsulation of [Dy(acac)₃(H₂O)₂] SMMs in multi-walled carbon nanotubes.³⁸ By modulating the diameter of the tubes, the encapsulation space can be fine-tuned to generate a desired ligand field around the Ln(III) ions, making them more promising for building novel hybrid devices. Rajaraman and co-workers have carried out a detailed theoretical investigation on encapsulating the [DyCp₂]⁺ family of complexes in SWCNTs, demonstrating that a high magnetic relaxation barrier can be retained upon encapsulation.³⁹

Computational tools are crucial for predicting mononuclear and polynuclear open-shell metal complexes’ geometry, electronic structure, and magnetic properties. Most importantly, computational tools are pivotal in designing and predicting next-generation SMMs.^6,40 Scalar relativistic density functional theory (SR-DFT) and multireference-based complete active space self-consistent field (CASSCF) methods,⁴¹ together with the treatment of dynamic correlations at the N-electron valence-state second-order perturbation theory (NEVPT2) level,⁴² are the main workhorse to compute the electronic structure and spin-Hamiltonian (SH) parameters of mononuclear transition metal and lanthanide complexes. In addition to calculating SH parameters, computational tools have been extensively employed to elucidate the underlying magnetic relaxation mechanism. Recently, the calculations of spin–phonon coupling have enabled us to gain a deeper understanding of the magnetic relaxation mechanism in mononuclear SIMs.^6,43–47 Periodic DFT calculations and ab initio molecular dynamics simulations have been used to shed light on the SMM behaviour in bulk.^48,49 Quantum chemical calculations have been widely applied to model complexes to explore their magnetic properties and understand their behaviour in hybrid assemblies. Some notable studies include the report of Chibotaru and co-workers, where a wide range of [LnX] and [LnX]⁺ model complexes (where Ln = Dy, Ho, Tb; X = O²⁻, S²⁻, Cl⁻, F⁻) were studied to explore the single-ion magnet behaviour.^50,51 Singh et al. studied linear [Ln(OH)_n] complexes (where Ln = Dy or Er) to investigate the influence of the geometry and symmetry around Ln(III) ions on magnetic anisotropy and relaxation.⁵⁴ Chibotaru and co-workers have thoroughly discussed the magnetic anisotropy associated with Ln(III) ions upon encapsulation in BNNTs and CNTs.^52,53 Considering the fact that the most robust axial limit can be achieved in simple linear [LnX] molecules (where X = non-metals),⁵¹ here we have explored the electronic structure, bonding, and magnetic anisotropy in rare-earth mononitrides [LnN] (where Ln = Dy and Tb), which is a promising neutral molecule that can offer giant magnetic anisotropy. Rare-earth mononitrides are refractory, durable ceramic materials utilized in various applications, including industrial coatings, catalysis, and superconductivity.⁵⁵ These [LnN] possess NaCl-type (Fm3m) structures and are noteworthy due to their unique semiconducting and ferromagnetic properties, making them relevant for spintronic devices.^56–58 Furthermore, a multitude of lanthanide nitride cluster fullerenes are reported, wherein three Ln(III) ions bonded to the nitride ion are encased within C₈₀ cages.^59,60 However, the isolation of molecular lanthanide mononitrides is scarce as they can directly interact with air and excess nitrogen.⁶¹ Nevertheless, findings indicate that molecular lanthanide nitrides have been captured in the gas phase through matrix-isolation infrared (IR) spectroscopy.^62,63

With very few reports on neutral isolated [LnN] molecules,^62,63 the present manuscript reports the study of the electronic structure, bonding, and magnetic anisotropy of [LnN] molecules (where Ln = Dy and Tb) and their encapsulation within 1D zigzag boron nitride nanotubes (BNNTs)⁶⁴ of varying diameters, specifically (8,0) and (9,0) to envisage novel hybrid assemblies with a giant magnetic blockade barrier. We have chosen boron nitride nanotubes (BNNTs) due to their several advantages over carbon nanotubes (CNTs), which include superior thermal and chemical stability, exceptional mechanical strength, high thermal conductivity and wide band gaps (∼5.5 eV) that are unaffected by chirality or morphology.^65,66 In the first section, we describe the bonding in [LnN] (where Ln = Dy(III) and Tb(III)) molecules using the CASSCF/NEVPT2 level of theory. In the second section, we describe the structural and energetic aspects of the perpendicular and parallel orientations of the [LnN] encapsulated in BNNT(8,0) and BNNT(9,0) rings using periodic density functional theory (pDFT). In the next section, we compute the magnetic anisotropy, g-factors, and crystal field parameters of naked [LnN] molecules and [LnN]@BNNT(8,0), [LnN]@BNNT(9,0) assemblies to assess their SMM behaviour. In the final section, we describe the effect of encapsulation on the magnetic anisotropy of [LnN] molecules and how these interactions could further tweak the SIM behaviour of the linear molecules.

Computational methodology

The cell and geometry optimizations of the pristine (8,0) and (9,0) nanotubes, [LnN], and [LnN]@BNNT(n,0) (where Ln = Dy(III), Tb(III) and n = 8, 9) were carried out using (pDFT) on the Vienna ab initio simulation package (VASP) 5.4.4 package.^67–69 We have used generalized gradient approximation (GGA) with the Perdew–Burke–Ernzerhof (PBE) exchange–correlation functional.⁷⁰ The plane wave energy cut-off E_cut was set to 600 eV for structural optimization. The Monkhorst–Pack grid 1 × 1 × 10 k-point has been used for the geometry optimization calculations, where the convergence threshold was set to 10⁻⁵ eV in energy and −0.05 eV A⁻¹ in force. Here, we used projector-augmented-wave (PAW) pseudopotentials for ionic core-electron interactions. Dispersion interactions were accounted for using Grimme's DFT D3 dispersion correction.⁷¹ The binding energy (B.E.) for the encapsulated system was computed using the following eqn (1):

B.E. = E_{LnN@BNNT(n,0)} − E_LnN − E_BNNT(n,0) … (1)

To further analyze the nature of interactions between the [LnN] and BNNT(8,0) and (9,0), we have performed PEDA using the BAND code.⁷² VASP-optimized geometries were exported to the BAND code for the PEDA and DOS calculations. All the PEDA calculations were performed with the PBE functional and 1 × 1 × 1 Γ-centered k-mesh, whereas a 1 × 1 × 3 regular k-point grid is used for the DOS calculations. Here, we employed Slater-type all-electron TZ2P basis sets^73,74 for Ln atoms and TZP for the rest of the atoms. Scalar relativistic effects were treated using the zeroth-order regular approximation (ZORA)⁷⁵ without any frozen core approximation. Dispersion corrections were applied using Grimme's D3 corrections with Becke–Johnson damping (D3BJ).⁷¹ We also performed extended transition state theory with natural orbitals for chemical valence (ETS-NOCV)⁷⁶ to further analyze the nature of the orbital interactions between the [LnN] and BNNTs. Here, we computed the interaction energy by fragmenting the hybrid system into two fragments: [LnN] as fragment 1 and BNNT(9,0)/(8,0) as fragment 2 (see the ESI‡ for details).

To compute the magnetic properties of the [LnN] molecule and [LnN]@BNNT(8,0)/(9,0) hybrid assemblies, we first prepared large cluster models. The cluster models were obtained by symmetrically cutting the BNNT fragment around the [LnN] molecule. Subsequently, the dangling bonds at the edges were saturated with the hydrogens, whose positions were optimized at the DFT level. Our model comprises three belts of B3N3 hexagons, with [LnN] sitting at the center. CASSCF calculations were performed on the cluster models to compute the magnetic anisotropy of [LnN] hybrid assemblies. All these calculations were carried out using the OpenMolcas code.⁷⁷ Here, we employed an active space of CAS(8,7) and computed 7 septets, 140 quintets, and 195 triplet spin-free states for the Tb(III) ion, while we computed 21 sextets for Dy(III) using an active space CAS(9,7). These 21 sextet states are good enough for the ⁶H_15/2 SOC spectrum and the magnetic properties of Dy(III)-based SIMs.^28,78–82 We have employed ANO-RCC type basis sets for all the atoms (see Table S7‡ for the contraction scheme). The Cholesky decomposition for two-electron integrals was adapted to improve storage efficiency. Finally, we have employed the SINGLE_ANISO module integrated into OpenMolcas to compute the g-values, crystal-field parameters, and wavefunction decomposition analysis. We have constructed the ab initio blockade barrier for magnetic relaxation analysis by computing the transition dipole moments between each KD.

To probe the nature of bonding in lanthanide nitrides, we have performed CASSCF calculations on the [DyN] analogue and developed the potential energy surface (PES) using ORCA/5.0.3.^83–85 Here, the active space is comprised of CAS(15,13), which includes nine electrons in seven 4f-orbitals, three Ln–N bonding orbitals (1σ + 2π), and the corresponding antibonding orbitals, i.e., six electrons in six orbitals. Upon converging CASSCF calculations, we carried out NEVPT2 employing the Douglas–Kroll–Hess (DKH2) approximation.⁸⁶ We have used the segmented all-electron relativistic contracted (SARC) triple-zeta basis set SARC-DKH2-TZVP⁸⁷ for the Ln(III) atom and the DKH-adapted def2-TZVP basis set for the N atom. We have also computed the IR spectra for the [LnN] molecules by carrying out optimization at the B3LYP level of theory with the SARC-DKH2-TZVP for Ln and DKH-def2-SVP for N. The dispersion interaction is addressed using Grimme's dispersion correction with the Becke–Johnson (D3BJ) scheme^71,88 implemented in ORCA. The potential energy scan is performed without any symmetry constraints.

Results and discussion

First, we developed the potential energy surface for the [DyN] molecule to shed light on the nature of bonding. Here, we have chosen an active space of CAS(15,13), which includes seven 4f-orbitals of Dy(III) and three bonding orbitals, i.e. one σ-bond (Dy(5d_z2)–N(2p_z)), two π-bonds (Dy(5d_xz))–N(2p_x) and (Dy(5d_yz)–N(2p_y)) and the corresponding antibonding orbitals. Using an active space of CAS(15,13), we scanned the Dy–N bond length from 1.70 Å to 1.90 Å with an increment of 0.02 Å (see Fig. 1a). The PES is constructed by computing two quartet and two doublet states in the CI procedure. Both the CASSCF and NEVPT2 calculations predicted that the minimum energy structure possesses a bond length of 1.764 Å, which matches closely with the previously reported DFT/CASSCF computed bond length of LnO/LnO⁺ molecules.^50,51 DKH-B3LYP optimization also yields a Dy–N bond distance of 1.750 Å, closely matching the NEVPT2 computed bond distance. Moreover, the computed stretching frequency of 824.5 cm⁻¹ at the DFT optimized geometry agrees with the previously isolated [DyN] in the argon matrix with the band at 810.5 cm⁻¹ (see Fig. S1‡).⁶² The spin-multiplicity sextet is observed to be the ground state throughout the scan. The ground state configuration can be written as (σ)²(π π′)⁴ (δ_4f)² (δ_4f)¹ (φ_4f)¹ (φ_4f)² (π_4f)¹ (π_4f)¹ (σ_4f)¹ (σ*)⁰ (π*π*′)⁰ (see Fig. 1b). A close inspection of the ground state electronic configuration yields the following population in the active orbitals, (σ)^1.93(π^1.93 π′^1.94) (δ_4f)^1.47 (δ_4f)^1.46 (φ_4f)^1.44 (φ_4f)^1.44 (π_4f)^1.10 (π_4f)^1.10 (σ_4f)^0.95 (σ*)^0.06 (π* ^0.06 π*′^0.06) (see Fig. 1b), suggesting a weakly correlated 4f-orbital in bonding. Moreover, the natural population analysis (NPA) on the NEVPT2 computed minima shows a gain of 0.90(0.87) in the vacant 6s(5d) orbitals of the Dy(III) ion, indicating their involvement in the bonding. Besides, we observed a loss of 0.13(1.83) from the 2s(2p) orbital of the N³⁻ ion, suggesting the involvement of these orbitals in the bonding.

Fig. 1 (a) Potential energy curves at the CASSCF and PT2 levels for the [DyN] bond length with the VASP-optimized bond lengths. (b) The bonding and antibonding orbitals are included in the active space (right) (isovalue 0.03 e bohr⁻³).

We observed a Mayer bond order of 2.69 for the lowest energy structure, indicating the presence of a triple bond between the Dy and N centers. The DFT/CASSCF computed [DyN] bond length is slightly smaller than the Dy–N bond distance observed from the pDFT calculations (1.856 Å).

Next, we encapsulated the [LnN] molecule in BNNTs of sizes (8,0) and (9,0) and performed the geometry optimization of the hybrid assemblies (see Computational details). Here, we have opted for two different orientations: (i) the parallel orientation, where the molecule lies along the axis of the nanotube, and (ii) the perpendicular orientation, where the [LnN] is oriented perpendicular to the diameter of the nanotube (see Fig. S2‡). We have named these assemblies as nLn_|| and nLn_⊥ (Ln = Dy(III) and Tb(III)), where [LnN] is encapsulated in BNNT(n,0) (n = 8, 9) in parallel and perpendicular fashion, respectively. Binding energy (B.E.) calculations predict the stabilization of these [LnN] molecules in the nanotube in parallel and perpendicular arrangements. The computed B.E. are ∼217.8(214.2) kJ mol⁻¹ for 8Dy_|| (8Tb_||), 183.9(180.2) kJ mol⁻¹ for 9Dy_|| (9Tb_||), 175.9(166.9) kJ mol⁻¹ for 8Dy_⊥ (8Tb_⊥) and 130.3(128.5) kJ mol⁻¹ for 9Dy_⊥ (9Tb_⊥), respectively (see Table S2‡). The perpendicular and parallel orientations of [LnN]@BNNT show sizable binding energy, with the parallel orientation being relatively more stable by ∼30–40 kJ mol⁻¹. Owing to the smaller diameter of BNNT(8,0) ∼6.4 Å compared to BNNT(9,0) ∼7.3 Å (see Table S1‡), we observed a strong interaction between the tube and [LnN] resulting in increased binding energy for [LnN]@BNNT(8,0). The computed [LnN] bond lengths are 1.876(1.880) and 1.878(1.879) Å for 8Dy_|| (8Tb_||) and 9Dy_|| (9Tb_||), respectively, and 1.901(1.906) and 1.857(1.858) Å for 8Dy_⊥ (8Tb_⊥) and 9Dy_⊥ (9Tb_⊥), respectively (see Fig. 2 and S3‡). A close inspection of structural parameters reveals that the B atom of the BNNT and the N atom of the [LnN] come much closer to each other ∼1.571 Å (exp. B–N bond length ∼1.4–1.6 Å) and form a weak bond in 8Dy_⊥ (8Tb_⊥), resulting in a giant binding energy for the perpendicular arrangement. Besides, the [LnN] molecule resides at the center of the tube in 8Dy_|| (8Tb_||) assemblies, whereas it moves closer toward one of the walls of the BNNT(9,0) in 9Dy_|| (9Tb_||) (see Fig. 2 and S3‡). Charge difference density (CDD) analysis shows a symmetric charge depletion around the [LnN] in 8Dy_|| (8Tb_||) assemblies, while an asymmetric depletion pattern is observed for 9Dy_|| (9Tb_||) (see Fig. 2). A comparison of non-covalent interaction (NCI)^89,90 plots of 8Dy_|| with 8Dy_⊥ shows the presence of a greenish region between the tube and [DyN] in 8Dy_||, indicative of stabilizing van der Waals interactions, while a reddish region is observed in 8Dy_⊥ which suggests a repulsive interaction between the tube and [DyN]. Moreover, the reddish region increases in 9Dy_|| compared to 8Dy_|| as the molecule moves closer to the wall of the BNNT(9,0) due to its larger diameter (see Fig. 2 and S3‡).

Fig. 2 pDFT-optimized structures of (a) 8Dy_||; (b) 8Dy_⊥; (c) 9Dy_|| and (d) 9Dy_⊥ with relevant bond parameters (bond distance in red represents the Dy–N bond distance of the [DyN] molecule); (e)–(h) CDD plots corresponding to optimized structures plotted with an isosurface = 0.00186 e Å⁻³; (i)–(l) non-covalent interaction (NCI) plots for the above structures plotted using Multiwfn⁹¹ and VMD software.⁹² Color code: Dy (green), N (blue), B (pink).

Next, we performed PEDA calculations to further analyze the interactions between the [LnN] and BNNTs. The calculated interaction energy (ΔE_int) consistently shows a stabilizing effect across all eight hybrid assemblies. Notably, the ΔE_int value is higher for [LnN]@BNNT(8,0) compared to [LnN]@BNNT(9,0), mirroring the trend seen in binding energy. As the complexes transition from 8Dy_|| (8Tb_||) to 9Dy_|| (9Tb_||), electrostatic and orbital interactions weaken due to the increased tube diameter, leading to a lower ΔE_int for [LnN]@BNNT(9,0). The 8Dy_⊥ and 8Tb_⊥ assemblies exhibit significantly high ΔE_int values of −427.3 kJ mol⁻¹ and −411.2 kJ mol⁻¹, nearly double those of the other assemblies. This remarkable stabilization is attributed to the formation of a B–N bond between the BNNT ring's B atom and the N atom in [LnN]. Further analysis of the orbital component of interaction energy (ΔE_orb) was conducted by decomposing it into NOCV donor/acceptor pairs (see the ESI‡ for details).

Next, we analyzed the density of states (DOS) of all the hybrid [LnN]@BNNT assemblies. The total density of states for all the [DyN]@BNNT assemblies are provided in Fig. 3, S6–S8.‡ The pristine BNNT shows no net spin at the Fermi level with a widely separated highest occupied crystal orbital (HOCO) and lowest unoccupied crystal orbital (LUCO), a typical manifestation of the wide band-gap semiconductors. For 9Dy_⊥, the projected DOS (pDOS) of the pristine [DyN] molecule shows a net magnetization (down-spin) at the Fermi level, which mainly arises from the unpaired electrons in the 4f orbitals of Dy(III), with minor contributions from the 4d and 5s orbitals near the Fermi level (see Fig. 3). Upon encapsulation of the [DyN] molecule in BNNT(9,0), we observed that 9Dy_⊥ displays net magnetization at the Fermi level, and the total DOS near the Fermi level closely resembles the pDOS of the [DyN] molecule. Moreover, the band gap in 9Dy_⊥ is smaller than that in the pristine BNNT(9,0), which indicates that the encapsulation helps enhance its conductivity by narrowing the band gap. A similar observation is made for the rest of the hybrid assemblies.

Fig. 3 Density of States plot for 9Dy_⊥. The green dotted line shows the Fermi level.

Magnetic anisotropy and single-ion magnet behaviour

Next, we analyzed the magnetic anisotropy of [LnN] and [LnN]@BNNT(8,0) and BNNT(9,0) assemblies in both the parallel and perpendicular orientations and compared them with pristine [LnN] units. All the calculations were performed on the cluster model by saturating the terminal dangling bonds with hydrogens. Here, we have studied the parallel and perpendicular systems separately in the BNNT(8,0) and BNNT(9,0) tubes and then made a comparison in the last section.

Magnetic anisotropy in [DyN] and [TbN]

CASSCF calculations were performed on the VASP-optimized geometry of [LnN] units where the Dy–N(Tb–N) bond distances are 1.856(1.855) Å, respectively(see Fig. 4). For [DyN], we observed the stabilization of a pure m_J |±15/2〉 with the eight low-lying KDs derived from the ⁶H_15/2 ground multiplet spanning over an energy range of 1821.4 cm⁻¹. The computed ground state g-values are completely Ising in nature (g_zz ∼ 20.011) with strictly no transverse components in the g-values. A purely axial environment in [DyN] maximizes the KD splitting with no mixing between the KDs, resulting in an Ising-type ground state. The principal g_zz axis of the g-tensor passes through the Dy–N bond. The first excited KD is purely m_J |±13/2〉 in nature, with strictly no mixing with other excited KDs (see Table S10‡). The computed transverse magnetic moment between the ground state KDs is extremely weak (6.44 × 10⁻²⁰μ_B), indicating a complete quenching of QTM within the ground state KDs. Analysis of the g-tensor orientation of all eight KDs shows strict collinearity, which indicates that magnetization reversal is likely to occur via KD8, resulting in the ab initio computed barrier height (U_cal) value of 1821.4 cm⁻¹ for [DyN]. The computed trends in the g-values, Ising-type ground state and barrier height of the [DyN] molecule match those of the previously studied [DyO]⁺ molecule.⁵¹

Fig. 4 Ab initio computed blockade barrier for (a) [DyN] and (b) [TbN] with the orientation of the g_zz axis (c). The yellow disk shows the beta-spin density.

For [TbN], we used an active space of CAS(8,7) and computed all 7 septets, 140 quintets, and 195 triplet spin-free states for the Tb(III) ion. SO-RASSI calculations reveal stabilization of the ⁷F₆ ground term with thirteen SOC states emerging from the ⁷F₆ ground term spanning over an energy range of 1820.9 cm⁻¹. For non-Kramers Tb(III) ions, the computed ground state g-values (g_zz = 18.014 with g_xx ∼ g_yy ∼ 0.000) are highly axial in nature. Due to the highly axial environment, we observed near degeneracy in the first twelve SOC states, arising from J = 6 to J = 1 with a negligible tunnel splitting (∼1 × 10⁻⁷ to 1 × 10⁻⁹ cm⁻¹), causing magnetization blockade up to the 12^th SOC state corresponding to the ⁷F₁ state. The magnetization reversal is likely to occur via the 13^th state, which is ⁷F₀ in nature and sets a record barrier height of ∼1821 cm⁻¹ in the [TbN] molecule. Interestingly, we observed that the energy span of the eight KDs in [DyN] and the 13 pseudo-KDs in [TbN] are of similar strength, highlighting that the complete axial ligand field generates nearly equal repulsion in the oblate-type f⁹ and f⁸ configurations. The [DyN] and the [TbN] molecules are highly anisotropic units exhibiting the record U_cal value of ∼1821 cm⁻¹, making them promising to explore for the encapsulation behavior.

Magnetic anisotropy in 8Dy_|| and 9Dy_||

Next, we performed calculations on 8Dy_|| (9Dy_||) to understand the magnetic anisotropy in the encapsulated systems and how the tube's diameter affects the single-ion anisotropy. In both the 8Dy_|| (9Dy_||) complexes, we observed stabilization of m_J |±15/2〉 as the ground state. The computed spin–orbit states corresponding to the ⁶H_15/2 ground multiplet span over an energy range of 1707.7 and 1666.8 cm⁻¹ for 8Dy_|| (9Dy_||), respectively (see Fig. 5 and S9‡). Compared to the [DyN] molecule, we observed that the energy span of the eight low-lying KDs is narrowed by 6.28% and 8.46% for 8Dy_|| (9Dy_||), respectively. These marginal changes in the SOC states arise due to the weak van der Waals interaction around the [DyN] molecule in both assemblies, as visualized in the NCI plots (see Fig. 2). The computed ground-state g-values are highly axial and display a typical Ising-type behavior (g_zz = 20.009, g_xx = g_yy = 0.000 for 8Dy_|| and g_zz = 20.007, g_xx = g_yy = 0.000 for 9Dy_||) for both the assemblies (see Table 1). Due to the highly axial environment in both 8Dy_|| and (9Dy_||) assemblies, we observed collinearity in the main magnetic axes up to the 7^th KD, suggesting that magnetization reversal is likely to occur via KD8, which is m_J |±1/2〉 in nature (see Tables S9 and S10‡). This suggests that the U_cal values are 1707 and 1666.8 cm⁻¹ for 8Dy_|| and 9Dy_||, respectively. In addition, we constructed a magnetic blockade barrier by computing the transverse magnetic moments between the KDs to analyze the mechanism of QTM. The computed ground state k_QTM value is 3.76 × 10⁻¹⁰μ_B for 8Dy_||, which is nearly ∼1000 times smaller than the 3.76 × 10⁻⁷μ_B observed for 9Dy_||. This k_QTM value difference is rooted in the changes in the ground-state transverse g-values in 8Dy_|| and 9Dy_||. The source of the increase in the transverse g-values and k_QTM is related to the change in the tube diameter. In 8Dy_||, we observed that the [DyN] molecules remain intact in the center of the tube and experience average electrostatic interactions from all sides of the wall. In the parallel orientation, the electrostatic interactions from the wall contribute to the equatorial ligand field. On the other hand, due to the large diameter of the BNNT(9,0) ring, the [DyN] molecule moves closer to one side of the wall of the tube (see Fig. 2c), which generates an unfavourable equatorial crystal field around the Dy(III) ion in 9Dy_|| resulting in a decrease in the KD span and an increase in the k_QTM values. To verify this, we have computed the crystal field parameters in 8Dy_|| and 9Dy_|| and compared them with those of the [DyN] molecule. In [DyN], we noticed that large axial B_k^q (q = 0) terms with negligible non-axial terms (q ≠ 0) support Ising-type g-values and negligible ground-state QTM. Besides, we observed dominant axial terms in 8Dy_|| and 9Dy_||; however, the values of the non-axial terms (q ≠ 0) significantly increased compared to the [DyN] molecule (see Table S12‡). A 10× decrease in the axial B_k^q (q = 0) parameters was observed for 9Dy_|| compared to 8Dy_|| because [DyN] experienced stronger interactions from non-axial directions when it moved closer to the wall side in the BNNT(9,0) rings. The decrease in the axial term of the CF arises due to the van der Waals interaction generated by the BNNT tubes in 8Dy_|| and 9Dy_|| perpendicular to the Dy–N bond length, which generates an equatorial ligand field. The decrease in barrier height is substantial for 9Dy_|| as the Dy–N molecule is much closer to the wall while marginal for 8Dy_||. A similar scenario was observed when neutral [DyO] was deposited on the hexagonal boron-nitride (h-BN) surface in a horizontal alignment, and the computed U_cal marginally changed compared to the pristine [DyO] molecule, owing to the large distance of 2.94 Å between the Dy(III) ion and the h-BN sheet.⁵⁰ In comparison to any other properties (g-values, spans of KDs), we observed that the crystal field parameters are susceptible to minor changes in the coordination environment which nicely explains why 9Dy_|| has a larger k_QTM value compared to 8Dy_||. The large barrier height in 8Dy_|| is further supported by the computed LoProp charges, which show marginally large charges on the nearest N atom of the tube in 9Dy_|| (−0.4156 a.u.) compared to 8Dy_|| (−0.4024 a.u.) (see Fig. S12‡). This slight increase in the LoProp charges indicates the stronger equatorial field offered upon encapsulation in 9Dy_||, decreasing the SOC states’ span and increasing the k_QTM values.

Fig. 5 (a) Orientation of the g_zz tensor for [DyN] encapsulated in BNNT(9,0) along with the blockade barrier for (a) parallel orientation and (b) perpendicular orientation. Color code: N (blue), B (pink), H (white), yellow (beta-spin density).

Table 1 The SINGLE_ANISO computed energies of the ground and excited Kramers doublets (KD), accompanied by their respective g-values, tunneling coefficients (k_QTM), and blocking barrier (U_cal)

	KD	g xx	g yy	g zz	k QTM	U cal (cm^−1)
DyN	KD1	0.000	0.000	20.011	6.44 × 10^−20	1821.4
	KD8	10.598	10.598	1.284	3.53 × 10
8Dy||	KD1	0.000	0.000	20.009	3.76 × 10^−10	1707.0
	KD8	10.713	10.496	1.293	3.53 × 10
8Dy⊥	KD41	0.000	0.000	20.005	1.26 × 10^−6	1669.5
	KD7	5.790	5.251	4.262	1.70 × 10
9Dy||	KD1	0.000	0.000	20.007	3.91 × 10^−7	1666.8
	KD8	1.013	5.319	15.261	3.58 × 10
9Dy⊥	KD1	0.000	0.000	20.008	1.27 × 10^−7	1939.7
	KD7	5.050	3.684	1.512	1.10 × 10

---

Magnetic anisotropy in 8Dy_⊥ and 9Dy_⊥

CASSCF calculations predict the stabilization of the m_J |±15/2〉 as the ground state in 8Dy_⊥ and 9Dy_⊥. The computed spin–orbital states corresponding to the ⁶H_15/2 ground multiplet span over an energy range of 1844.7(1989.1) cm⁻¹ for 8Dy_⊥ (9Dy_⊥), respectively (see Table S8‡). Compared to the [DyN] molecule, the span of the eight low-lying KDs is expanded by ∼1.27% (9.21%) in the case of 8Dy_⊥ (9Dy_⊥). The calculated ground state g-values are highly axial and display a typical Ising type behavior (g_zz = 20.009, g_xx = g_yy = 0.000 (8Dy_⊥) and g_zz = 20.008, g_xx = g_yy = 0.000 (9Dy_⊥)) for both the assemblies. The first excited KD is purely m_J |±13/2〉 in nature and located ∼533(569) cm⁻¹ higher in energy from the ground state KD for 8Dy_⊥ (9Dy_⊥), respectively (see Table 1). The principal z-axis of the g-tensor in 8Dy_⊥ and 9Dy_⊥ lies along the Dy–N bond vector. A careful analysis of the g_zz orientation of all eight low-lying KDs suggests that the g_zz orientation is collinear up to KD5 for 8Dy_⊥ and KD6 for 9Dy_⊥. Analysis of the g-tensor orientation hints that the magnetic relaxation will occur from KD6 in 8Dy_⊥ and KD7 in 9Dy_⊥, leading to the U_cal values of 1669.5 and 1939.7 cm⁻¹ for 8Dy_⊥ and 9Dy_⊥, respectively (see Fig. 5b). The computed transverse magnetic moments between the ground state KDs are 1.26 × 10⁻⁶ (1.27 × 10⁻⁶)μ_B for 8Dy_⊥ (9Dy_⊥), indicative of a complete blockade of ground-state QTM. The computed U_cal indicates that 9Dy_⊥ shows a higher barrier height than 8Dy_⊥ and [DyN] molecules. This is in line with previous theoretical studies indicating an increase in the U_cal value for the [DyO]⁺ molecule deposited perpendicular to the MgO surfaces due to the enhancement in the axial crystal field effect emerging from the negatively charged O atom at a shorter distance between [DyO]⁺ and the MgO surface.⁵¹ The smaller diameter of the (8,0) tube is unfavourable for hosting the [DyN] molecule in perpendicular mode as it results in the formation of a covalent bond between the [DyN] molecule and the B atom of the tube, resulting in a decreased SIM performance for 8Dy_⊥ compared to 9Dy_⊥, due to the slight increase in the Dy–N bond length in 8Dy_⊥ (1.901 Å) compared to 9Dy_⊥ (1.856 Å). Moreover, the B–N bond formation results in charge transfer from the N atom of the [DyN] molecule to the coordinated B atom, which led to much smaller LoProp charges on the N atom in 8Dy_⊥ than in 9Dy_⊥. The elongated Dy–N bond length coupled with weaker LoProp charges on the N atom of [DyN] decreases the axiality in 8Dy_⊥, which directly affects the KD span and g-values (see Fig. S12‡). Based on the results obtained, we observed the following trend in U_cal values: 9Dy_⊥ > 8Dy_|| > 8Dy_⊥ ∼ 9Dy_||. Most importantly, all the [DyN]@BNNT(n,0) (n = 8, 9) assemblies show robust SIM performance, irrespective of the orientation. The trends obtained in U_cal clearly suggest that size and orientation are equally important in determining the performance of Dy(III)-based single-ion magnets when encapsulated in the nanotube. Moreover, it is important to note that further increasing the ring size beyond BNNT(9,0) for the perpendicular orientation does not guarantee an increase in U_cal values as the crystal field effects weaken with the increase in distance. This has nicely been observed for the [DyO]⁺@C₆₀, where the magnetic blockade barrier marginally changes when the Dy(III) atom moves away from the hexagonal ring of C₆₀.⁵¹

Magnetic anisotropy in 8Tb_|| and 9Tb_||

Next, we performed calculations on 8Tb_|| and 9Tb_|| to compute the SOC states, SH parameters, and tunnel splitting values to shed light on the magnetic anisotropy of [TbN] upon encapsulation (see Table S9‡). In both 8Tb_|| and 9Tb_|| assemblies, we find the ⁷F₆ as the ground state with thirteen m_J states spanned over ∼1626.5 and ∼1535.7 cm⁻¹ for 8Tb_|| and 9Tb_||, respectively (see Table S9‡). Compared to the [TbN] molecule, we observed that the splitting of the thirteen m_J states corresponding to the ⁷F₆ term is narrowed by ∼10% and ∼15% for 8Tb_|| and 9Tb_||, respectively. The extent of decrease in the SOC splitting upon encapsulation of the [TbN] in 8Tb_|| and 9Tb_|| is much higher than the [DyN] molecule. The computed ground state g-values are highly axial (g_xx = g_yy = 0.000 and g_zz = 17.939 for 8Tb_|| and g_xx = g_yy = 0.000 and g_zz = 17.938 for 9Tb_||) with stabilization of m_J |±6〉 in both the assemblies (see Table 2, and S9‡). In both cases, the computed g_zz orientation is aligned towards the Tb–N bond vector. For 8Tb_||, we observed a ground state tunnel splitting of ∼1 × 10⁻⁸ cm⁻¹, which is extremely small and predicts negligible ground state QTM for a non-Kramers Tb(III) ion. The computed ground-state tunnel splitting for 8Tb_|| matches closely the [TbN] value, which indicates that parallel encapsulation results in minor changes in the ground state QTM values. For 8Tb_||, we observed that the tunnel splitting is minimal up to KD3 m_J |±4〉 (<1 × 10⁻⁴ cm⁻¹), which indicates that the magnetic relaxation is likely to occur via the 4^th pseudo-KD with a tunnel splitting value of 6.73 × 10⁻² cm⁻¹. This predicts a U_cal value of ∼1015 cm⁻¹ for 8Tb_||. Surprisingly, a giant tunnel splitting of 1.37 × 10⁻³ cm⁻¹ is observed for 9Tb_||, indicating the presence of the dominant QTM within the ground state pseudo-KD. As observed in 9Dy_||, the [TbN] also exhibits closer proximity to the wall of the tube, which generates an unfavourable crystal field and consequently increases the tunnel splitting. The computed CF parameters show that the axial B_k^q (q = 0) terms are negative and of nearly similar strength in [TbN], 8Tb_|| and 9Tb_||, indicating axiality in all the cases. In contrast, the computed non-axial B_k^q (q ≠ 0) CF terms have grown significantly in 8Tb_|| and 9Tb_|| compared to the [TbN] value, indicating the presence of QTM in both assemblies (see Table S13‡). Analysis of the Δ_tun values of ground and excited pseudo-KDs suggests that magnetic relaxation is likely to occur via the pseudo-KD3, resulting in a U_cal value of ∼674.7 cm⁻¹ for 9Tb_||, which is nearly ∼33% of the U_cal value observed for 8Tb_||. In a similar manner to 8Dy_||, we observed that the [TbN] molecule encapsulated in the BNNT(8,0) helps in maintaining the higher axiality, which leads to a giant blockade barrier for magnetic relaxation. In a complete axial ligand field environment, we observed that the span of low-lying KDs/pseudo-KDs is the same for the [TbN] and [DyN] molecules. For 8Dy_|| and 9Dy_||, we observed a change in the U_cal value by <10% compared to the [DyN] molecule. Owing to similar structural changes in 8Tb_|| and 9Tb_||, here, we observed that the U_cal value changed by 30–65% compared to the [TbN] molecule, which highlights that the magnetic anisotropy and SIM behaviour in [TbN] molecule changes drastically reduced upon encapsulation.

Table 2 The SINGLE_ANISO-computed energies of the ground and excited Kramers doublets (KDs), accompanied by tunnelling splitting (k_QTM), by their respective g-values, blocking barrier (U_cal)

	Pseudo KD	Δ tun	g xx	g yy	g zz	U cal (cm^−1)
TbN	KD1	0.00 × 10	0.000	0.000	18.014	1820.9
8Tb||	KD1	7.10 × 10^−9	0.000	0.000	18.011	1015.8
	KD4	6.73 × 10^−2	0.000	0.000	8.071
8Tb⊥	KD1	8.05 × 10^−5	0.000	0.000	17.937	350.5
	KD2	1.09 × 10^−2	0.000	0.000	14.555
9Tb||	KD1	1.37 × 10^−3	0.000	0.000	18.010	674.7
	KD3	2.42 × 10^−1	0.000	0.000	11.269
9Tb⊥	KD1	6.22 × 10^−5	0.000	0.000	17.941	808.9
	KD3	1.62 × 10^−2	0.000	0.000	11.122

---

Magnetic anisotropy in 8Tb_⊥ and 9Tb_⊥

For 8Tb_⊥ (9Tb_⊥) assemblies, we find that the CASSCF computed energies of the ⁷F₆ ground term span over an energy range of ∼1821.0(1848.3) cm⁻¹, which is much larger than that for 8Tb_|| and 9Tb_|| (see Table S9‡). Compared to the [TbN] molecule, the computed span of low-lying pseudo-KDs is expanded by ∼1.5% for 9Tb_⊥, while it remains the same for 8Tb_⊥. This indicates that the ligand field axiality in [TbN] is complemented upon encapsulation in the BNNT(9,0). The computed ground state g-values are highly axial in nature (g_xx = g_yy = 0.000 and g_zz = 17.937 for 8Tb_⊥ and g_xx = g_yy = 0.000 and g_zz = 17.941 for 9Tb_⊥) with stabilization of the m_J |±6〉 in both the assemblies (see Table 2 and S9‡). The computed tunnel splitting for the ground state pseudo-KD m_J |±6〉 is found to be extremely small for both 8Tb_⊥ (9Tb_⊥), with a value of 8.05 × 10⁻⁵(6.22 × 10⁻⁵) cm⁻¹. Compared to the 9Tb_||, we observed a ∼1000 times smaller k_QTM value for the 9Tb_⊥. This large difference in the observed k_QTM value arises from the structural differences where the tube acts as an axial ligand to the Tb(III) ion in the perpendicular arrangement. On the other hand, we observed an opposite trend in BNNT(8,0), where the k_QTM value for the 8Tb_⊥ is ∼1000 times higher than the k_QTM value for 8Tb_||. The perpendicular arrangement of the Tb(III) ion in BNNT(8,0) led to the formation of a bond between the N atom of the [TbN] and the B atom of the tube, which led to an increase in the Tb–N bond length and decreases the axiality in Tb(III) ion compared to 9Tb_⊥. Due to the high Δ_tun value in the ground and first excited pseudo-KDs, magnetic relaxation will likely occur via the first excited pseudo-KD in 8Tb_⊥, resulting in a U_cal value of ∼350.5 cm⁻¹. On the other hand, magnetic relaxation occurs via the third pseudo-KD in 9Tb_⊥, resulting in a U_cal value of ∼809 cm⁻¹.

After careful analysis of the magnetic anisotropy and magnetic relaxation of [LnN] in BNNT(8,0) and BNNT(9,0), it is evident that the diameter of the BNNT(8,0) suits well for the encapsulation of parallel [LnN] and minimal structural changes result in a giant magnetic anisotropy for 8Dy_||/8Tb_||. In 8Dy_||, we observed the U_cal value of 1707 cm⁻¹, which is only ∼7% reduced compared to the [DyN] molecule. A similar scenario is observed for 8Tb_||; however, the extent of the decrease is relatively higher than that of the Dy analogue. In addition, we observed that the large diameter of the BNNT(9,0) easily accommodates the [LnN] molecule in a perpendicular arrangement without any drastic structural changes. The axiality of the [LnN] molecule is further leveraged by the surrounding atoms of the BNNT(9,0), which led to a sizable increase of ∼100 cm⁻¹ in 9Dy_⊥ compared to the [DyN] molecule.

Conclusion

In summary, this study offers a comprehensive exploration of the structure, bonding, and magnetic anisotropy of highly anisotropic [LnN] (Ln = Tb(III) and Dy(III)) molecules, as well as their encapsulation within zigzag boron nitride nanotubes of sizes BNNT(8,0) and BNNT(9,0) to shed light on the magnetic anisotropy of hybrid assemblies relevant for the development of new generation magnetic materials. NEVPT2 calculations predict the presence of a triple bond in [Ln-N], where the bonding mainly originated from one σ-bond (Dy(5d_z2)–N(2p_z)) and two π-bonds Dy(5d_xz)–N(2p_x) and Dy(5d_yz)–N(2p_y), with the 4f orbital being weakly correlated and not contributing towards the bonding. Our in-depth analysis of structural and magnetic properties suggests that lanthanide mononitride possesses a giant barrier height >1800 cm⁻¹. Upon encapsulation, we observed that the alignment of the [LnN] molecules parallel to the nanotube's axis is energetically more favourable in BNNT(8,0) tubes as the smaller diameter of the BNNT(8,0) tube generates dominant dispersion interactions. For the Kramers-type Dy(III) ion, we observed that all four studied hybrid assemblies are qualified for the best SIMs. However, BNNT(9,0) outperforms the rest in the series owing to the highly axial ligand field provided by the BN atoms of the tube in perpendicular mode, rendering a high barrier height of 1939 cm⁻¹. On the other hand, the tunnel splitting for non-Kramers Tb(III) ions is highly sensitive towards the crystal field environment, and only 8Tb_|| is qualified to be a good SIM. The energetic preference and single-ion anisotropy analysis suggests that both 8Dy_|| and 8Tb_|| are outstanding SIMs, thanks to the smaller diameter of the BNNT(8,0) tube, which keeps the molecule intact in the centre. Our findings suggest that BNNTs are excellent hosts for encapsulating highly anisotropic [LnN] molecules that can display a giant magnetic relaxation barrier. In conclusion, we described the importance of encapsulation behaviour and the possibility of tuning the tube diameter and molecular orientation to maximize the SIM behaviour while integrating Ln(III)-based SIMs to generate strongly anisotropic hybrid assemblies.

Data availability

The data supporting this article have been included as part of the ESI.‡

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

SKS acknowledges the Science and Engineering Research Board (CRG/2023/002936) and IIT Hyderabad for their generous funding. KK and SM acknowledge PMRF. SKS acknowledges the PARAM Seva Computing Facility under the National Supercomputing Mission at IIT Hyderabad.

References

1. D. N. Woodruff, R. E. P. Winpenny and R. A. Layfield, Chem. Rev., 2013, 113, 5110–5148.
2. R. E. P. Winpenny, Angew. Chem., Int. Ed., 2008, 47, 7992–7994.
3. L. Bogani and W. Wernsdorfer, Nat. Mater., 2008, 7, 179–186.
4. L. Ungur and L. F. Chibotaru, Chem. – Eur. J., 2017, 23, 3708–3718.
5. J. D. Rinehart and J. R. Long, Chem. Sci., 2011, 2, 2078.
6. I. Tarannum, S. Moorthy and S. K. Singh, Dalton Trans., 2023, 52, 15576–15589.
7. S. K. Singh, C. J. Cramer and L. Gagliardi, Inorg. Chem., 2020, 59, 6815–6825.
8. S. Moorthy, I. Tarannum, K. Kumari and S. K. Singh, Dalton Trans., 2024, 53, 12073–12079.
9. K. Kumari and S. K. Singh, Dalton Trans., 2024, 53, 16495–16511.
10. F.-S. Guo, B. M. Day, Y.-C. Chen, M.-L. Tong, A. Mansikkamäki and R. A. Layfield, Science, 2018, 362, 1400–1403.
11. C. A. P. Goodwin, F. Ortu, D. Reta, N. F. Chilton and D. P. Mills, Nature, 2017, 548, 439–442.
12. A. H. Vincent, Y. L. Whyatt, N. F. Chilton and J. R. Long, J. Am. Chem. Soc., 2023, 145, 1572–1579.
13. P. Evans, D. Reta, G. F. S. Whitehead, N. F. Chilton and D. P. Mills, J. Am. Chem. Soc., 2019, 141, 19935–19940.
14. J. Liu, Y.-C. Chen, J.-L. Liu, V. Vieru, L. Ungur, J.-H. Jia, L. F. Chibotaru, Y. Lan, W. Wernsdorfer, S. Gao, X.-M. Chen and M.-L. Tong, J. Am. Chem. Soc., 2016, 138, 5441–5450.
15. X.-L. Ding, Y.-Q. Zhai, T. Han, W.-P. Chen, Y.-S. Ding and Y.-Z. Zheng, Chem. – Eur. J., 2021, 27, 2623–2627.
16. Y.-C. Chen, J.-L. Liu, L. Ungur, J. Liu, Q.-W. Li, L.-F. Wang, Z.-P. Ni, L. F. Chibotaru, X.-M. Chen and M.-L. Tong, J. Am. Chem. Soc., 2016, 138, 2829–2837.
17. Y. Chen, J. Liu, Y. Lan, Z. Zhong, A. Mansikkamäki, L. Ungur, Q. Li, J. Jia, L. F. Chibotaru, J. Han, W. Wernsdorfer, X. Chen and M. Tong, Chem. – Eur. J., 2017, 23, 5708–5715.
18. Y. Chen, J. Liu, W. Wernsdorfer, D. Liu, L. F. Chibotaru, X. Chen and M. Tong, Angew. Chem., Int. Ed., 2017, 56, 4996–5000.
19. A. B. Canaj, S. Dey, E. R. Martí, C. Wilson, G. Rajaraman and M. Murrie, Angew. Chem., Int. Ed., 2019, 58, 14146–14151.
20. S. Jia, X. Zhu, B. Yin, Y. Dong, A. Sun and D. Li, Cryst. Growth Des., 2023, 23, 6967–6973.
21. S. K. Gupta, T. Rajeshkumar, G. Rajaraman and R. Murugavel, Chem. Sci., 2016, 7, 5181–5191.
22. C. Zhao, Z. Zhu, X.-L. Li and J. Tang, Inorg. Chem. Front., 2022, 9, 4049–4055.
23. M. Li, J. Han, H. Wu, Y. Du, Y. Liu, Y. Chen and S. Chen, Dalton Trans., 2022, 51, 16964–16972.
24. J. Acharya, N. Ahmed, J. F. Gonzalez, P. Kumar, O. Cador, S. K. Singh, F. Pointillart and V. Chandrasekhar, Dalton Trans., 2020, 49, 13110–13122.
25. P. Kalita, N. Ahmed, S. Moorthy, V. Béreau, A. K. Bar, P. Kumar, P. Nayak, J.-P. Sutter, S. K. Singh and V. Chandrasekhar, Dalton Trans., 2023, 52, 2804–2815.
26. A. B. Canaj, S. Dey, C. Wilson, O. Céspedes, G. Rajaraman and M. Murrie, Chem. Commun., 2020, 56, 12037–12040.
27. S. Bala, G.-Z. Huang, Z.-Y. Ruan, S.-G. Wu, Y. Liu, L.-F. Wang, J.-L. Liu and M.-L. Tong, Chem. Commun., 2019, 55, 9939–9942.
28. P. Kalita, K. Kumari, P. Kumar, V. Kumar, S. K. Singh, G. Rogez and V. Chandrasekhar, Dalton Trans., 2024, 53, 10521–10535.
29. W.-J. Xu, Q.-C. Luo, Z.-H. Li, Y.-Q. Zhai and Y.-Z. Zheng, Adv. Sci., 2024, 11, 2308548.
30. G. Serrano, L. Poggini, M. Briganti, A. L. Sorrentino, G. Cucinotta, L. Malavolti, B. Cortigiani, E. Otero, P. Sainctavit, S. Loth, F. Parenti, A.-L. Barra, A. Vindigni, A. Cornia, F. Totti, M. Mannini and R. Sessoli, Nat. Mater., 2020, 19, 546–551.
31. V. Corradini, A. Candini, D. Klar, R. Biagi, V. Renzi, A. Rizzini, N. Cavani, U. Pennino, S. Klyatskaya, M. Ruben, E. Velez-Fort, K. Kummer, N. Brookes, P. Gargiani, H. Wende and M. Affronte, Nanoscale, 2018, 10, 277–283.
32. K. Katoh, J. Sato, R. Nakanishi, F. Ara, T. Komeda, Y. Kuwahara, T. Saito, B. K. Breedlove and M. Yamashita, J. Mater. Chem. C, 2021, 9, 10697–10704.
33. M. Del Carmen Giménez-López, F. Moro, A. La Torre, C. J. Gómez-García, P. D. Brown, J. Van Slageren and A. N. Khlobystov, Nat. Commun., 2011, 2, 407.
34. V. Tangoulis, N. Lalioti, J. Parthenios, N. Boukos, O. Malina, J. Tuček and R. Zbořil, Inorg. Chem., 2018, 57, 6391–6400.
35. R. Nakanishi, J. Satoh, K. Katoh, H. Zhang, B. K. Breedlove, M. Nishijima, Y. Nakanishi, H. Omachi, H. Shinohara and M. Yamashita, J. Am. Chem. Soc., 2018, 140, 10955–10959.
36. A. A. Popov, S. Yang and L. Dunsch, Chem. Rev., 2013, 113, 5989–6113.
37. A. A. Popov and L. Dunsch, J. Am. Chem. Soc., 2007, 129, 11835–11849.
38. R. Nakanishi, M. A. Yatoo, K. Katoh, B. K. Breedlove and M. Yamashita, Materials, 2017, 10, 7.
39. R. Nabi, R. K. Tiwari and G. Rajaraman, Chem. Commun., 2021, 57, 11350–11353.
40. P. Kumar, J. F. Gonzalez, P. P. Sahu, N. Ahmed, J. Acharya, V. Kumar, O. Cador, F. Pointillart, S. K. Singh and V. Chandrasekhar, Inorg. Chem. Front., 2022, 9, 5072–5092.
41. B. O. Roos, P. R. Taylor and P. E. M. Sigbahn, Chem. Phys., 1980, 48, 157–173.
42. S. Guo, M. A. Watson, W. Hu, Q. Sun and G. K.-L. Chan, J. Chem. Theory Comput., 2016, 12, 1583–1591.
43. A. Mattioni, J. K. Staab, W. J. A. Blackmore, D. Reta, J. Iles-Smith, A. Nazir and N. F. Chilton, Nat. Commun., 2024, 15, 485.
44. D. Aravena, J. Phys. Chem. Lett., 2018, 9, 5327–5333.
45. A. Lunghi, F. Totti, R. Sessoli and S. Sanvito, Nat. Commun., 2017, 8, 14620.
46. A. Lunghi, F. Totti, S. Sanvito and R. Sessoli, Chem. Sci., 2017, 8, 6051–6059.
47. A. Lunghi and S. Sanvito, Sci. Adv., 2019, 5, eaax7163.
48. A. Lunghi, M. Iannuzzi, R. Sessoli and F. Totti, J. Mater. Chem. C, 2015, 3, 7294–7304.
49. N. F. Chilton, Inorg. Chem., 2015, 54, 2097–2099.
50. W. Zhang, A. Muhtadi, N. Iwahara, L. Ungur and L. F. Chibotaru, Angew. Chem., Int. Ed., 2020, 59, 12720–12724.
51. L. Ungur and L. F. Chibotaru, Inorg. Chem., 2016, 55, 10043–10056.
52. D. Liu, X. Guo, X. Zhang, A. A. Al-Kahtani and L. F. Chibotaru, Inorg. Chem., 2024, 63, 3769–3780.
53. D. Liu, X. Guo, X. Zhang, A. A. Al-Kahtani and L. F. Chibotaru, Dalton Trans., 2023, 52, 11243–11253.
54. S. K. Singh, T. Gupta and G. Rajaraman, Inorg. Chem., 2014, 53, 10835–10845.
55. A. R. H. Preston, S. Granville, D. H. Housden, B. Ludbrook, B. J. Ruck, H. J. Trodahl, A. Bittar, G. V. M. Williams, J. E. Downes, A. DeMasi, Y. Zhang, K. E. Smith and W. R. L. Lambrecht, Phys. Rev. B:Condens. Matter Mater. Phys., 2007, 76, 245120.
56. C.-G. Duan, R. F. Sabirianov, W. N. Mei, P. A. Dowben, S. S. Jaswal and E. Y. Tsymbal, J. Phys.: Condens. Matter, 2007, 19, 315220.
57. J. P. Evans, G. A. Stewart, J. M. Cadogan, W. D. Hutchison, E. E. Mitchell and J. E. Downes, Phys. Rev. B, 2017, 95, 054431.
58. K. Kneisel, M. Maddah, J. Chan, Y. Xu, C. Casey-Stevens, K. Van Koughnet, W. Holmes-Hewett, H. J. Trodahl and F. Natali, ACS Omega, 2024, 9(48), 47842–47847.
59. S. Stevenson, J. P. Phillips, J. E. Reid, M. M. Olmstead, S. P. Rath and A. L. Balch, Chem. Commun., 2004, 2814–2815.
60. T. Zuo, C. M. Beavers, J. C. Duchamp, A. Campbell, H. C. Dorn, M. M. Olmstead and A. L. Balch, J. Am. Chem. Soc., 2007, 129, 2035–2043.
61. M. Takano, A. Itoh, M. Akabori and K. Minato, J. Phys. Chem. Solids, 2005, 66, 697–700.
62. S. P. Willson and L. Andrews, J. Phys. Chem. A, 1999, 103, 1311–1321.
63. G. V. Chertihin, W. D. Bare and L. Andrews, J. Phys. Chem. A, 1998, 102, 3697–3704.
64. D. Golberg, Y. Bando, C. C. Tang and C. Y. Zhi, Adv. Mater., 2007, 19, 2413–2432.
65. X. Blase, A. Rubio, S. G. Louie and M. L. Cohen, EPL, 1994, 28, 335.
66. C. Zhi, Y. Bando, T. Terao, C. Tang, H. Kuwahara and D. Golberg, Adv. Funct. Mater., 2009, 19, 1857–1862.
67. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
68. G. Kresse and J. Furthmüller, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
69. G. Kresse and D. Joubert, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
70. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
71. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
72. G. te Velde and E. J. Baerends, Phys. Rev. B:Condens. Matter Mater. Phys., 1991, 44, 7888–7903.
73. E. Van Lenthe and E. J. Baerends, J. Comput. Chem., 2003, 24, 1142–1156.
74. M. A. Watson, N. C. Handy and A. J. Cohen, J. Chem. Phys., 2003, 119, 6475–6481.
75. E. Van Lenthe, E. J. Baerends and J. G. Snijders, J. Chem. Phys., 1994, 101, 9783–9792.
76. M. P. Mitoraj, M. Parafiniuk, M. Srebro, M. Handzlik, A. Buczek and A. Michalak, J. Mol. Model, 2011, 17, 2337–2352.
77. F. Aquilante, J. Autschbach, A. Baiardi, S. Battaglia, V. A. Borin, L. F. Chibotaru, I. Conti, L. De Vico, M. Delcey, I. F. Galván, N. Ferré, L. Freitag, M. Garavelli, X. Gong, S. Knecht, E. D. Larsson, R. Lindh, M. Lundberg, P. Å. Malmqvist, A. Nenov, J. Norell, M. Odelius, M. Olivucci, T. B. Pedersen, L. Pedraza-González, Q. M. Phung, K. Pierloot, M. Reiher, I. Schapiro, J. Segarra-Martí, F. Segatta, L. Seijo, S. Sen, D.-C. Sergentu, C. J. Stein, L. Ungur, M. Vacher, A. Valentini and V. Veryazov, J. Chem. Phys., 2020, 152, 214117.
78. K. Bernot, J. Luzon, L. Bogani, M. Etienne, C. Sangregorio, M. Shanmugam, A. Caneschi, R. Sessoli and D. Gatteschi, J. Am. Chem. Soc., 2009, 131, 5573–5579.
79. G. Cucinotta, M. Perfetti, J. Luzon, M. Etienne, P.-E. Car, A. Caneschi, G. Calvez, K. Bernot and R. Sessoli, Angew. Chem., Int. Ed., 2012, 51, 1606–1610.
80. M. Briganti, G. F. Garcia, J. Jung, R. Sessoli, B. L. Guennic and F. Totti, Chem. Sci., 2019, 10, 7233–7245.
81. M. S. Norre, C. Gao, S. Dey, S. K. Gupta, A. Borah, R. Murugavel, G. Rajaraman and J. Overgaard, Inorg. Chem., 2020, 59, 717–729.
82. D. Shao, W.-J. Tang, Z. Ruan, X. Yang, L. Shi, X.-Q. Wei, Z. Tian, K. Kumari and S. K. Singh, Inorg. Chem. Front., 2022, 9, 6147–6157.
83. F. Neese, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2012, 2, 73–78.
84. F. Neese, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2022, 12, e1606.
85. F. Neese, F. Wennmohs, U. Becker and C. Riplinger, J. Chem. Phys., 2020, 152, 224108.
86. T. Nakajima and K. Hirao, Chem. Rev., 2012, 112, 385–402.
87. D. A. Pantazis and F. Neese, J. Chem. Theory Comput., 2009, 5, 2229–2238.
88. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
89. E. R. Johnson, S. Keinan, P. Mori-Sánchez, J. Contreras-García, A. J. Cohen and W. Yang, J. Am. Chem. Soc., 2010, 132, 6498–6506.
90. A. Otero-de-la-Roza, E. R. Johnson and J. Contreras-García, Phys. Chem. Chem. Phys., 2012, 14, 12165–12172.
91. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
92. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.

---

Footnotes

† This work is dedicated to Professor V. Chandrasekhar on the occasion of his 65^th birthday.

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

---