Cubic aromaticity in ligand-stabilized doped Au superatoms

The magnetic response of valence electrons in doped gold-based [ M @Au 8 L 8 ] q superatoms (M = Pd, Pt, Ag, Au, Cd, Hg, Ir, and Rh; L = PPh 3 ; and q = 0, + 1, + 2) is studied by calculating the gauge including magnetically induced currents (GIMIC) in the framework of the auxiliary density functional theory. The studied systems include 24 different combinations of the dopant, total cluster charge, and cluster structure (cubic-like or oblate). The magnetically induced currents (both diatropic and paratropic) are shown to be sensitive to the atomic structure of clusters, the number of superatomic electrons, and the chemical nature of the dopant metal. Among the cubic-like structures, the strongest aromaticity is observed in Pd-and Pt-doped [ M @Au 8 L 8 ] 0 clusters. Interestingly, Pd-and Pt-doping increases the aromaticity as compared to a similar all-gold eight-electron system [ Au 9 L 8 ] + 1 . With the recent implementation of the GIMIC in the deMon2k code, we investigated the aromaticity in the cubic and butterfly-like M@Au 8 core structures, doped with a single M atom from periods 5 and 6 of groups IX–XII. Surprisingly, the doping with Pd and Pt in the cubic structure increases the aromaticity compared to the pure Au case not only near the central atom but encompassing the whole metallic core, following the aromatic trend Pd > Pt > Au. These doped (Pd, Pt)@Au 8 nanoclusters show a closed shell 1S 2 1P 6 superatom electronic structure corresponding to the cubic aromaticity rule 6 n + 2.


I. INTRODUCTION
The ubiquitous concept of aromaticity 1 has been popularized since the discovery of benzene by Faraday in 1825 2 and Kekulé's proposal for its structure in 1865. 3Even though aromaticity has been widely used to explain the stability in molecules and coordination compounds, 4,5 related with spectroscopy 6 and magnetic properties, 7 and employed to describe chemical reactivity 8,9 in several compounds, there is no consensus 10 about its nature due to the lack of a quantum mechanical observable to directly measure it.
The high symmetry of organic molecules or compounds is commonly related to aromaticity, as is the stability of their electronic structure related to the closed-shell structure. 11,12Although there have been several efforts to establish a criterion to determine aromaticity in the different bare and protected metal clusters [13][14][15][16][17] and successful attempts have been made to indirectly correlate its effects with measurable quantities such as chemical shifts in nuclear magnetic resonance (NMR; 18 e.g., the aromatic nature of a compound can be observed through the down-field shift in the 1 H probe NMR caused by magnetic shielding 2,7 ), there are still some questions to address.Indeed, the cubic geometry has been explored in different systems to find a remarkable stability by relating the electron delocalization and the electronic structure with a concept named cubic aromaticity.To mention some examples, there have been studies about this cubic aromaticity in Zn 8 , 19 Mn 8 , 20 Co 13 O 8 , 21 and CBe 8 H 12 22 clusters.
Aromaticity has also been studied in doped monolayerprotected gold superatoms. 23These systems are of interest because of the promising applications that can enhance targeted physicochemical properties, such as luminescence and catalysis, 24,25 by the The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp][28] Consequently, the effects of the closed shell electronic structure and the composition of the metallic core on the physico-chemical properties have been explored.0][31] The preferred sites were observed to be at the center location for Pd and Pt and on the surface of the core for Cd and Hg. 32The doping effect has also been studied for the elements of the group IX, Rh, and Ir. 33Similarly, as for Pd and Pt, the endohedral location was preferred for the dopants in M@Au 24 (SR) clusters.The hydrogen evolution reaction was studied in doped M@Au 24 nanoclusters where M = Pd and Pt. 34Doping was found to decrease the adsorption free energy for hydrogen, thus improving the catalytic activity of the doped clusters, with Pt doped clusters having the best performance.
An additional example of a widely studied system are the superatoms with a metal core based on Au 9 .Of these compositions, clusters without doping protected by PPh 3 have been crystallized in three different symmetries: butterfly-like with D 2h symmetry, 35 crown-like with approximately D 4d symmetry, 36,37 and metastable cubic with O h symmetry recently obtained by the crystallization of a metastable state. 38As a consequence of the interest in tuning the composition of the Au metal cores, clusters M@Au 8 with Pd 39 or Pt 40 at the center have also been experimentally observed.These clusters can be used in the synthesis of larger nanoclusters via hydride doping 41 and have possible applications in catalysis, 42 motivating studies of their formation 43 and further modification. 44n this work, we have systematically studied the magnetic response of valence electrons of 24 different doped gold-based [M@Au 8 L 8 ] q superatoms (M = Pd, Pt, Ag, Au, Cd, Hg, Ir, and Rh; L = PPh 3 ; and q = 0, +1, +2) having 6-10 delocalized superatom electrons in the metal core, with the aim of using the magnetically induced diatropic current as the descriptor for aromaticity. 45We have employed our recent implementation of the gauge-including magnetically induced currents (GIMIC) 46 in the framework of the auxiliary density functional theory (ADFT).The implementation uses the gauge-independent atomic orbitals (GIAO) introduced by London 47 to solve the gauge origin problem in the computation of the magnetic properties of molecules and takes advantage of efficient recurrence relations to compute the shielding tensor, 48 resulting in a high performance of the computation of the magnetically induced current susceptibility.We show that the magnetically induced current is sensitive to the atomic structure of the clusters, the number of superatomic electrons, and the chemical nature of the dopant metal.

A. Response of electrons to an external magnetic field
When a diamagnetic object (molecule and nanocluster) is under the influence of an external magnetic field ⃗ B, the movement of the electrons generates an induced current that is a sum of diamagnetic and paramagnetic contributions.The induced current can flow in two directions relative to the external magnetic field: a classical diatropic one (clockwise) and a nonclassical paratropic one (anticlockwise).Molecules or nanostructures that experience a strong diatropic response are called aromatic, those that sustain mainly paratropic currents are called antiaromatic, and those that have almost the same strengths for both contributions are called non-aromatic. 49hus, a formally called magnetic aromaticity/antiaromaticity must be related to the ability of a ring or arrangements of atoms in a closed chain to sustain a global induced current (clockwise/counterclockwise direction) in the presence of an external perpendicular magnetic field.This induced global current circulation is commonly known as the ring-current model. 50,51ensity functional theory (DFT) computations have been used to indicate the molecular response properties under the influence of an external magnetic field (e.g., on the magnetic susceptibilities 52 and chemical shifts 48 ).In fact, these properties are computed as second order energy change (ΔE (2) ) in a molecule/nanoparticle interacting with the external magnetic field.Thus, the magnetic shielding tensor σ C,ηλ (η, λ = x, y, z) at the position of the nucleus C can be computed as the second derivative of the energy with respect to the external magnetic field ⃗ B and the nuclear magnetic moment Another second order energy expression involves a divergenceless vector potential ⃗ where ⃗ J(⃗ r ) is the electron current density.Then, substituting the explicit form of the vector potential 46 and differentiating with respect to the Cartesian components of the vectors ⃗ μC and ⃗ B, we get the nuclear magnetic shielding tensor via the Biot-Savart law, 50,[54][55][56][57][58][59][60] where εηατ(η, α, τ = x, y, z) is the Levi-Cività tensor and ⃗ C is the position of the nucleus C (within the implementation, we employ a.u., and then, μ 0 /4π = 1/c 2 ).The quantity J (λ) τ is defined as It is an element of the so-called current density tensor (CDT), 59,61,62 which characterizes the response of the current density to an external magnetic field.λ indicates differentiation with respect to the corresponding Cartesian coordinate of the magnetic field.Finally, the magnetically induced current (MIC) density can be calculated as the first order correction to the total current density.The MIC only includes direct dependence on basis functions, first derivatives of basis functions, and density matrix and magnetically perturbed density matrix.In this work, those quantities are computed in a standard SCF procedure in deMon2k. 63he Journal of Chemical Physics ARTICLE scitation.org/journal/jcp

B. Computational details
Ground state DFT calculations were performed starting from crystallographic data (cube 38 and butterfly 35 ) using the real-space grid based GPAW program. 64,65The frozen core approximation was used for all the elements, with 15 valence electrons for Rh and Ir, 16 for Pd and Pt, 17 for Ag, 11 for Au, 12 for Cd, and 18 for Hg.The projector augmented wave (PAW) setups for heavy metals included scalar-relativistic corrections.The structures were optimized employing the generalized gradient approximation (PBE) exchange-correlation functional. 66The optimization was finished when the residual forces on each atom were less than 0.05 eV/Å.The real-space grid spacing was 0.2 Å.To identify the symmetries of the superatomic orbitals, the Y lm analysis 67 was performed for each cluster around the center of mass with a cut-off radius of 4.0 Å.
The magnetically induced current density was computed by the recent implementation 46 in deMon2k 63 employing gauge-including magnetically induced currents theory.The Stuttgart-Dresden (SDD) pseudopotential 68 (Rh and Ir-17 valence electrons, Pd and Pt-18 valence electrons, Ag and Au-19 valence electrons, and Hg and Cd-20 valence electrons) along with the DZVP 69 basis set (H, C, and P).PBE 66 was used to consider the exchange and correlation effects.All the computations were performed in combination with the GEN-A2 * auxiliary function set.

III. RESULTS AND DISCUSSION
The clusters studied in this work with composition [M@Au 8 (PPh 3 ) 8 ] q (q = 0, +2: M = Cd, Pd, Pt, and Hg; and q = +1: M = Rh, Ag, Ir, and Au) were obtained by doping the corresponding all-gold clusters: the nearly perfect body centered cubic (BCC) structure 38 and the approximate D 2h symmetry derived from an icosahedron, 35 as depicted in Fig. 1.These structures were chosen for their close similarity.Considering only the metal core, the transformation from the cubic shape to the butterfly shape can be achieved by compressing the distances of the Au atoms located on the yz plane in the y-coordinate and elongating them in the z-coordinate.For the sake of simplicity, we employ a notation for the two different structures of clusters [M@Au 8 L 8 ] q=0,+1,+2 (L = PPh 3 ) as M@1 q=0,+1,+2 for the cubic structures and M@2 q=0,+1,+2 for the butterfly-like structures, as shown in Fig. 1.In both cases, the Au based structures were doped with the single atom (Rh, Ir, Pd, Pt, Ag, Cd, and Hg) occupying the central position.All the hypothetical systems were optimized, as discussed in Sec.II B. The total energies, energy differences, and HOMO-LUMO (H-L) gaps are summarized in Tables S1 and S2.

A. Geometrical description
The distances (r in Å) from the central atom to the eight surrounding Au atoms can be easily seen in Fig. 2. All the M@1 q structures preserved an approximately cubic structure after the geometry optimization, where small distortions from the perfect BCC structure were observed.In the same way, all the M@2 q butterfly-like structures preserved the D 2h symmetry, with the main differences with respect to the pure Au 9 core located in the longest distance from the central atom.To verify the different distortions with respect to the pure Au 9 cases, the principal moments of inertia were computed for all the structures M@1 q and M@2 q , as summarized in Table S3, and the ratio with respect to the highest moment of inertia I 1 in Table S4.As shown in the bottom panel of Fig. 2, the shapes of the metallic core are not perfect cubes in any doping structure neither in the pure Au 9 case where the distances range from 2.68 to 2.70 Å (see Table S5).However, in any case, the structures were distorted, and they were still not far away from the BCC structure, as demonstrated by the ratio of the moments of inertia (Table S4).For the cubic structures, the shortest distances were observed for Pd@1 0 (2.64-2.66Å) with Hg@1 +2 being the structure that has the most expanded structure (2.72-2.74Å) (Table S5).In the case of the butterfly-like structures, two sets (4:4) of characteristic distances were shown (Table S5).The shortest ones were observed in Rh@2 +1 and the longest in Hg@2 +2 .

B. Electronic structure
A criterion for the closed electronic structures is determined by the number of valence electrons (n) confined in the superatomic core given by the following equation first proposed by some of us: 67 FIG. 2. Distances from the central atom to the surrounding eight Au atoms for the cubic (bottom row) M@1 q=0,+1,+2 and butterfly-like (top row) M@2 q=0,+1,+2 structures.The doping element and charge are indicated in different colors.

ARTICLE scitation.org/journal/jcp TABLE I.
Closed electronic structure configuration for the lowest minimum energy of M@(1, 2) q .Valence electrons provided by the dopant atom M (ν), total charge (q), and number of superatom electrons (n).
where νx is the number of valence electrons provided by the xth constituent metal atom, N is the number of metal atoms, j is the number of one-electron-withdrawing ligands, and q is the total charge.The number of valence electrons ν donated by the different dopant atoms and the superatomic electron count according to Eq. ( 5) can be seen in Table I.The superatomic electron configuration is also shown for the most stable isomer (1 or 2) of each cluster.These configurations were based on the Y lm analysis, 67 which was performed for the clusters to distinguish the spherical symmetries of the frontier orbitals, as visualized in Figs.S1-S4.In the following, we discuss the electronic structures based on the chemical characterization (ν) from the dopant atom.

Ir/Rh @ Au 8
As shown in Table I, Ir and Rh have ν = −1, making them electron-withdrawing dopants that diminish the superatomic electron count.Thus, the doped clusters with a charge of +1 are sixelectron superatoms.Due to the splitting of the P orbitals, all these clusters have relatively large H-L gaps (0.89-1.14 eV).The (Ir, Rh)@2 +1 isomers have large H-L gaps and lower energies than their cubic counterparts.For the (Ir, Rh)@2 +1 clusters, HOMO and HOMO-1 are not superatomic orbitals, but localized d orbitals of the central dopant atom.LUMO in both, however, is the 1Px orbital upshifted due to the Jahn-Teller effect.In the (Ir, Rh)@1 +1 isomers, the superatomic electron configuration is 1S 2 1P 4 .However, some of the frontier orbitals on the occupied side show d-type symmetry and have a strong contribution from the orbitals from the dopant.

Pd/Pt @ Au 8
Pd and Pt provide zero valence electrons to the superatom, making them neutral dopants.The (Pd,Pt)@1 0 clusters are thus eight-electron superatoms, with the shell-closing 1S 2 1P 6 and significant H-L gaps of 1.43 eV and 1.51 eV, respectively.In the (Pd, Pt)@2 0 cluster, the HOMO and HOMO-1 are localized d orbitals of the central atom, as shown in Fig. S2.The H-L gap is reduced compared to the cubic isomers, which are also lower in energy.
In the (Pd, Pt)@2 +2 clusters with the 1S 2 1P 4 superatomic electron configuration, the Jahn-Teller effect causes the 1P orbitals to split into three non-degenerate orbitals.The 1Pz orbital (see Fig. S2) is downshifted due to the elongation of the cluster in the z-direction (lowest moment of inertia I 3 ; see Table S3) and the 1Px is upshifted (highest moment of inertia I 1 ; see Table S3).This results in the large H-L gap of these six-electron superatoms, as shown in Fig. 3(b).The cubic isomers (Pd,Pt)@1 +2 experience somewhat smaller splitting due to the more symmetric shape of the metal core and consequently also have a smaller H-L gap.The cubic isomers are also less stable, as shown in Fig. 3(a).

Ag/Au @ Au 8
Au and Ag both have ν = +1, making the doped clusters of q = +1 eight-electron superatoms.(Au, Ag)@1 +1 both have a configuration of 1S 2 1P 6 .For the Ag cluster, the three P orbitals are almost degenerate, whereas for Au, 1Pz is a little higher in energy, reflecting the slightly larger distribution in the principle moments of inertia, as shown in Table S4.The electronic structure of (Au, Ag)@2 +1 can be written as 1S 2 1P 4 y,z 1D 2 z 2 .Here again, the 1Px orbital is significantly upshifted and also the 1D 2 z orbital downshifted due to the shape of the metal core, causing the change in the filling order of the superatomic orbitals.As a consequence, the H-L gap remains large,

ARTICLE
scitation.org/journal/jcpalthough smaller than that for the cubic isomers.The lower-energy isomer is 1 for Au and 2 for Ag.
4. Cd/Hg @ Au 8 Since Cd and Hg are donors of valence electrons with ν = +2, the clusters with q = 0 have a superatom electron count of n = 10.(Cd, Hg)@1 0 has practically no H-L gap, and in fact, the D-like frontier orbitals have fractional occupation numbers.These isomers are higher in energy than (Cd, Hg)@2 0 , which have 1S 2 1P 4 y,z 1D 2 z 2 1P 2 x configurations and exhibit small H-L gaps.The 1Px and 1D 2 z orbitals are shifted similarly as for (Au, Ag)@2 +1 .

C. Magnetically induced current densities
To study the magnetically induced current density ⃗ J (1) (⃗ r) in all the systems, an external magnetic field ⃗ B = (0, 0, 1T) was imposed in the z-direction (Fig. 1) for all the clusters.To visualize the currents in different sectional planes of the clusters, the current densities were projected on three x-y planes at z = 0.0 Å, z = 0.75 Å, and z = 2.22 Å, starting from the z-coordinate of the central atom.Additionally, the paratropic, diatropic, and total contributions to these projected current densities were circularly integrated at different radii r from the center of the plane (0, 0, z), showing which contribution dominates at each radius.The projected current density fields and the circularly integrated currents for all the systems are shown in Figs.S5-S7.
To quantify the total current contributions over each of the three selected planes, radial integrals of the circularly integrated currents were calculated.These contributions are summarized in Tables S6-S8.We want to focus on those systems that showed the lowest energy, as summarized in Table I.
For (Rh, Ir)@2 +1 , the major diatropic contribution is localized in the z = 0.0 Å plane (Fig. S5) with the total radial integral contributions of −2.79 and −2.22 Å 2 nA/T (Table S6), respectively.Then, in these butterfly-like systems, the central atom is behaving as a single atom exhibiting a strong localized ring current within ∼1.0 Å from the nucleus.The z = 2.22 Å plane shows also a localized diatropic contribution near the positions of the Au atoms (Fig. S5) that is accompanied by a paratropic contribution within r = 1.0 ∼ 2.0 Å from the center of the plane.In both systems, a global induced current is not observed.
A similar behavior is observed for the (Cd, Hg)@2 0,+2 systems where a strong localized ring current is present near the central atom, higher for the neutral systems than for the charged ones (Fig. S7), with a total radial integral contribution of −1.71 and −2.34 Å 2 nA/T for Cd and Hg, respectively (see Table S8).
Finally, for the three other systems that showed a minimum in the butterfly-like structure, (Pd, Pt)@2 +2 and Ag@2 +1 , a localized diatropic contribution is again observed around the central doping atom (Figs.S6 and S5) with total radial integral contributions of −1.94 and −1.93 Å 2 nA/T for the Pd and Pt clusters, respectively (Table S7), and −1.48 Å 2 nA/T for the Ag cluster (Table S6).
The only systems having the cubic structure lower in energy than the butterfly-like one are Pd@1 0 , Pt@1 0 , and Au@1 +1 , which FIG. 4. Vectorial field of the magnetically induced current density and the circularly integrated current density for the cubic eight-electron superatoms Pd@1 0 , Pt@1 0 , and Au@1 +1 at the z-coordinate (Å).Diatropic (blue) and paratropic (red) contributions are shown.The total integrated current is shown in black.The external magnetic field B is pointing toward the reader.-S7; in Å 2 nA/T) at radius r = 4.0 Å from the center of an x-y plane located at (0, 0, 0) in clusters M@(1) 0,+1,+2 .
Comparing the cubic (1) and butterfly (2) structures in Figs.S5-S7, we can see that while both support diatropic currents at the z = 0 plane, this contribution is localized only around the central atom in butterfly structures.However, the cubic isomers have a substantial diatropic current over the whole metal core, as shown in Table II.Tables S5-S7 show that the total radially integrated current contributions of the cubic clusters are consistently more diatropic than in their butterfly-like counterparts.In both structures, the total integrated diatropic current also increases with the superatom electron count.This seems natural since a larger delocalized electron count can create a stronger response current in the same volume.Interestingly, the radial integrals for cubic structures show also diatropic contributions on the other planes, indicating a global induced current that is not diatropic only at the center but extends to the whole metallic core.This indicates that structures 1 are aromatic.The superatom electron count is shown in Table II Figure 4 shows the effect of Pd and Pt doping for the induced currents compared to the cubic all-Au cluster.Both the maximum circularly integrated diatropic current and the total radially integrated diatropic contribution are bigger for the Pd and Pt doped clusters than for the all-Au cluster (for instance, in the plane of the dopant, the total diatropic contributions are as follows: Pd: −6.05, Pt: −5.59, and Au: −3.09 Å 2 nA/T).This fact indicates that these dopants increase the aromaticity of the cluster.The order of aromaticity of these structures induced by the central atom follows the trend Pd > Pt > Au.

IV. CONCLUSIONS
In this work, we have studied the electronic structure and magnetically induced currents of doped [M@Au 8 L 8 ] q superatoms (M = Pd, Pt, Ag, Au, Cd, Hg, Ir, and Rh; L = PPh 3 ; and q = 0, +1, +2) having 6-10 delocalized superatom electrons in the metal core.We employed the gauge-including magnetically induced currents framework to compute the current susceptibility tensor and then the magnetically induced current density.This scheme has been previously used to give a quantitative measurement of the aromaticity/antiaromaticity in molecules via diatropic and paratropic contributions. 46Analysis of all 24 different systems showed the following trends: (i) total integrated diatropic currents increase with the superatom electron count, (ii) the cubic-like clusters have generally larger diatropic currents than the butterfly-like clusters, and (iii) the total integrated diatropic currents, hence aromaticity, increases for eightelectron cubic clusters as M = Au < Pt < Pd, i.e., doping with zerovalent metals Pd and Pt increases the aromaticity.We hope that this work encourages further studies about the correlation between stability in closed shell superatoms and magnetic response properties, such as the magnetically induced current, helping to reach a better understanding about the aromatic character of complex 3D metal clusters and nanostructures.

SUPPLEMENTARY MATERIAL
See the supplementary material for total energies, HOMO-LUMO energy gaps, principal moments of inertia, PDOS, and MICs of all the systems studied in this work.

FIG. 3 .
FIG. 3. (a) The relative energies for the two isomers for all the different dopings and charge states.The numbers and the lengths of the vertical bars indicate the difference in energy between the isomers.The vertical position of the bars has no physical meaning.(b) The HOMO-LUMO gaps for the same systems.

TABLE II .
Total radial integrals of the circularly integrated current strength (Figs.S5