Predictive optical photoabsorption of Ag₂₄Au(DMBT)₁₈⁻ via efficient TDDFT simulations

Atomically precise nanoclusters, also called nanomolecules or monolayer-protected clusters (MPCs), i.e., metal nanoclusters protected by a layer of coating ligands with a well-defined stoichiometry and chemical structure, constitute a class of materials of great interest in fundamental science and applications^1–27 and are also ideal systems to test and validate our physicochemical understanding of the phenomena exhibited by metal nanostructures at work.^28–31 Knowledge of the MPC precise chemical entity allows researchers to apply multiple characterization schemes, both experimental and theoretical, to these systems, thus cross-validating methods and results in a rigorous framework. This possibility is especially important in the field of the computational prediction of the structure and properties of materials. Given the availability of widely different computational methods and schemes employing different approximations and techniques, it is crucial to assess precisely the performances and also the pitfalls of these different choices, so as to identify the computational protocols best reproducing experimental results and properties. These protocols can then be employed to obtain accurate information also in those cases in which theory is called to provide missing information that would be difficult or even impossible to obtain experimentally. In other words, once validation is successfully achieved, the computational protocol can be trusted within well-understood limits and can be applied with grounded confidence in a predictive way to design new materials with desired optical properties.

The experimental photo-absorption spectrum of MPC compounds is one of the properties in strong need of validation studies.^32,33 Although qualitative agreement between experimentally observed and theoretically simulated spectra has been established,^34–36 quantitative agreement has been achieved only in few cases.^37–39 For example, to the best of our knowledge, the first response calculations using a hybrid exchange-correlation (xc-)functional both in the Kohn–Sham part and in the time-dependent density-functional theory (TDDFT) kernel were reported in Ref. 37, where real-time TDDFT/B3LYP simulations were conducted on the “green gold” Au₃₀(S-tBu)₁₈ nanomolecule. In that work, via a detailed comparison with low-temperature experimental spectra, it was clearly demonstrated that a sophisticated and computationally expensive level of theory (hybrid xc-functionals in both Kohn–Sham orbitals and TDDFT kernel and experimentally validated geometry) is mandatory to predict, with good accuracy, both the position of the low-energy peak and especially the ratio of oscillator strengths between the low-energy peak and the higher energy excitations of Au₃₀(S-tBu)₁₈, where both the position of the low-energy peak and the intensity ratio give this aliphatic MPC its unique green color.

The challenge in achieving quantitative accuracy in theoretical predictions is due to the difficulties brought about by the co-existence in these species of three components with very different electronic properties: (1) the metal core where the metallic bond is dominant, (2) the cationic-like metal species in the “staple” motifs interacting simultaneously with the metal core and the ligands, and (3) the anionic-like dichalcogenide atoms in the ligands, again simultaneously interacting with the metal core surface and the cationic-like metal species. The simultaneous presence of these different M–M^δ+, X^δ−–M^δ+, and X^δ−–M interactions with partial ionic and partial covalent character (where M stands for a metal atom and X for a dichalcogenide atom) and the sensitivity of the property of interest (here, the photo-absorption spectrum) to fine details of the atomistic structure make the choice of a theoretical approach able to predict accurately both the structure and the chiroptical response of these compounds very difficult.^37,38 For example, it has been shown that, at the time-dependent density-functional theory (TDDFT) level, simulations using non-hybrid xc-functionals, such as gradient-corrected ones (GGA), in their kernel tend to underestimate the absorption threshold of MPCs.³⁴ This in turn poses serious questions on the reliability or at least the accuracy of these approaches in the prediction of more complex, photo-absorption-derived phenomena, such as optical photoluminescence, in which the low-energy excitations close to the absorption threshold play an essential role.^40–43 

To make progress in this field, here we focus on the photo-absorption spectrum of the Ag₂₄Au(DMBT)₁₈ singly negative anion, where DMBT is the 2,4-dimethylbenzenethiolate ligand. This is a complex and realistic case exhibiting many of the prototypical features and physicochemical phenomena of the MPC class of materials. Indeed, the Ag₂₄Au(DMBT)₁₈⁻ compound has the characteristics of presenting: (1) both Ag and Au atoms, (2) aromatic thiolate ligands directly attached to the sulfur atoms, and (3) a well-defined and robust structure. These characteristics are advantageous for a validation study because they allow one to test: (1) an alloyed system, (2) ligands electronically conjugated with the metal core, and (3) an unequivocal configurational framework. The present study thus significantly extends our previous study on the Ag₂₅(CH₃)₁₈⁻ anion, in which we considered an MPC exhibiting a single metal element (Au) and aliphatic (non-aromatic and non-sulfur-conjugated) ligands³⁸ and our previous work on the aliphatic Au₃₀(S-tBu)₁₈ nanomolecule.³⁷ To study the more complex Ag₂₄Au(DMBT)₁₈⁻ compound, we apply advanced, recently developed simulation and analysis tools^44,45 that allow us to predict the chiro-optical spectrum of this medium-sized cluster at the TDDFT level using a hybrid xc-functional in the kernel, to analyze it in terms of single-particle excitations, and to compare the results with those obtained using a non-hybrid xc-functional. Additionally, we also compare and validate theoretical predictions with experimental photoabsorption spectra obtained in the range from room to low (77 K) temperatures. Clearly, MPCs with conjugated ligands, as considered here, are more challenging than MPCs with aliphatic ligands, as previous literature regarding the accuracy of density-functional theory (DFT) and TDDFT in the prediction of geometries⁴⁶ and optical spectra⁴⁷ of conjugated organic molecules has shown that a high-level of theory is needed for these components, and thus, in turn, it is also necessary to make predictions for MPCs with conjugated ligands in which the conjugated organic ligands resonating with the metal–sulfur core must be described accurately to achieve a comparable level of accuracy for the composite system. However, we stress that qualitatively analogous conclusions hold also for MPCs with aliphatic ligands,^37,38 as we discuss at the end of Sec. IV.

As an outcome of this study, (i) we first demonstrate that the combination of hybrid xc-functionals in the kernel and experimental geometry is validated as a quantitatively predictive tool and that we can (ii) identify the reasons for the failure of non-hybrid TDDFT in achieving accurate simulations of the optical response of MPC compounds and (iii) single out which type of excitation is not correctly modeled by the use of non-hybrid approaches, thus providing insights into the nature of electronic excitations in these materials. The so-derived information should be useful in future studies, both computational and experimental, as it provides progress in methodological tools and physical understanding.

Density-functional theory (DFT) is nowadays the most popular approach in quantum chemistry (QC) due to the optimal balance of accuracy and computational effort.⁴⁸ Although in principle an exact theory, in practice its accuracy depends significantly on the approximation of the xc-functional employed in any given case and is also system-dependent.⁴⁹ The most basic and also the most computationally efficient of these approximations is the Local Density Approximation (LDA).⁴⁸ At a higher level, but still typically requiring a comparable computational effort, lie the gradient-corrected xc-functionals, among which the generalized gradient approximation, GGA, and in particular the Perdew–Burke–Ernzerhof (PBE) flavor, is the most used.⁵⁰ A substantial increase in computational complexity is instead entailed by so-called hybrid xc-functionals, i.e., xc-functionals including a component of the Hartree–Fock exchange, among which B3LYP is one among the most used variants.^51,52 Since the Hartree–Fock exchange is non-local, the evaluation of its matrix elements within hybrid xc-functionals represents a significant challenge for all QC codes based on periodic basis sets, such as plane waves, as well as localized but non-transferable basis sets, such as Slater-type orbitals (STOs). A compromise between GGA and hybrid xc-functionals is provided by mixed schemes, such as the LB94 one,⁵³ in which a correction term exhibiting an asymptotic Coulombic tail is added to the KS Hamiltonian and corrects one (however, not all) of the issues of semi-local xc-functionals.

DFT is a ground-state theory. Its response analogue is the Time-Dependent DFT (TDDFT) approach.⁵⁴ The most common implementation of the TDDFT approach consists in expanding the molecular Kohn–Sham (KS) orbitals as linear combinations of atomic functions and recast the TDDFT problem in terms of diagonalizing a matrix omega (⁠$Ω$⁠) along the density-matrix formulation of Casida.⁵⁵ In practice, the KS equations are solved first,

Then, the $Ω$ matrix is diagonalized,

where the matrix elements of the $Ω$ matrix can be calculated from the KS orbitals and eigenvalues (1),

and in Eq. (3) the coupling matrix elements K_ia,jb are defined as

In Eq. (4), $fστXC$ is the XC kernel, with σ and τ being spin variables. The Casida problem consists in finding eigenvalues and eigenvectors of the matrix $Ω$ [Eq. (2)]. The eigenvalues of Eq. (2) correspond to squared excitation energies, while from the eigenvectors the intensities (absolute oscillator strengths) can be extracted.⁵⁵ Note, importantly, that the choice of the xc-functional translates within TDDFT into a double choice of: (1) which xc-functional to use to derive KS orbitals and (2) which xc-functional to use in the XC kernel $fστXC$ in Eq. (4).

Even using efficient diagonalization algorithms, such as Davidson’s, it becomes very difficult to calculate photoabsorption spectra over a wide excitation energy range when big systems are considered. The Davidson algorithm is efficient on large $Ω$ matrices, but only for the extraction of a relatively small number of lowest eigenvalues and eigenvectors. Therefore, the Casida TDDFT algorithm cannot be employed in practice to calculate photoabsorption spectra over a wide energy range. To overcome this problem, we have proposed the complex polarizability TDDFT (polTDDFT) algorithm.⁵⁶ This approach is able to treat very large systems and avoid the bottleneck of the diagonalization via a direct solution of the response equations and provides very good accuracy when contrasted with the exact Casida solution.^38,57 The reader is referred to the original work for a detailed description of the algorithm⁵⁶ and its implementation in the ADF program.⁵⁸ 

In practice, the photoabsorption spectrum $σω$ is calculated point by point, from the imaginary part of the dynamical polarizability $αω$⁠,

This expression is of practical interest when the polarizability is calculated for complex frequency, i.e., ω = ω_r + iω_i, where the real part ω_r is the scanned photon frequency (energy) and ω_i is the imaginary part, which corresponds to a broadening of the discrete lines and can be interpreted as a pragmatic inclusion of the excited state finite lifetime. The complex dynamical polarizability is calculated by solving the following non-homogeneous linear system, working with a basis set of density-fitting functions,

In Eq. (6), S is the overlap matrix between fitting functions, b is the unknown vector with the expansion coefficients $bμω$ of $ρz(1)$⁠, and d is the frequency-dependent vector corresponding to the known non-homogeneous term, and finally, the elements of the frequency-dependent matrix M are

In Eq. (7), χ_KS refers to the Kohn–Sham frequency-dependent dielectric function and K to the kernel. Note that the matrix element in Eq. (7) is between density fitting functions, and therefore, the present implementation of the algorithm allows one to employ only density-dependent kernels. For this reason, the density-matrix-dependent Hartree–Fock exchange kernel, which is an ingredient of hybrid xc-functionals, cannot directly be employed in polTDDFT. However, it is still possible to employ hybrid kernels with non-local Hartree–Fock exchange components in polTTDDFT, if the Hybrid Diagonal Approximation (HDA) is employed. The HDA that we have recently proposed (see Ref. 44 to which we refer for further details) is a very efficient and robust approximation to TDDFT: by correcting only the diagonal matrix elements of the Casida $Ω$ matrix in Eq. (2) with non-local exchange terms, while including local-density (Coulomb-type) terms of the kernel in the off-diagonal elements, the computational complexity is drastically reduced while keeping an excellent accuracy.⁴⁴ The HDA has been implemented in both Casida and polTDDFT methods within the AMS suite of codes [https://www.scm.com/product/ams/] and allows users an important saving in the computational effort as well as, crucially, the possibility of using hybrid xc-functionals, such as B3LYP in polTDDFT, thereby achieving an accuracy otherwise impossible using LDA or GGA xc-funtionals, as we will report below.

The cluster geometry has been taken from the experiment,⁵⁹ then it has been fully or partially optimized in the staple units (both cases have been considered), where in the partial optimizations, the metal core and the staple units (i.e., the metal and sulfur atoms) have been kept fixed at their experimental geometry, whereas only the metal core has been frozen in the fully relaxed staple geometries. Geometry relaxations were performed using the CP2K package⁶⁰ at the DFT/PBE level⁵⁰ with the addition of the Grimme-D3 dispersion terms.⁶¹ DFT simulations, based on the hybrid Gaussian/Plane-Wave scheme (GPW),⁶² employed Double-Zeta-Valence-plus Polarization (DZVP) basis sets⁶³ to represent the DFT Kohn–Sham orbitals, GTH pseudopotentials for describing the core-electrons of all the atomic species,⁶⁴ and an auxiliary set of plane-waves whose cutoff was set to 300 Ry.

Once geometries were thus derived as described above, optical response calculations were performed with the Amsterdam Modeling Suite (AMS) [https://www.scm.com/product/ams/] set of programs,^65,66 employing a basis set of Slater-type orbital (STO) functions of Triple Zeta plus Polarization (TZP) quality. Scalar relativistic effects have been included via the Zero Order Regular Approximation⁶⁷ (ZORA). Three different xc-functionals have been considered: (1) the standard GGA PBE,⁵⁰ (2) LB94 that carries the correct asymptotic Coulombic tail,⁵³ which is an important feature to obtain more accurate TDDFT results but is not supported by standard GGA functionals, and (3) the hybrid B3LYP that includes a portion of Hartree–Fock non-local exchange and provides consistently the most accurate predictions. As recalled in Sec. II, the use of hybrid xc-functionals is computational demanding because of the non-local nature of the Hartree–Fock exchange especially in the TDDFT kernel, which can be efficiently managed when Gaussian basis functions are used, but becomes dramatically inefficient when using an STO basis. We can, however, include hybrid xc-functionals via the Hybrid Diagonal Approximation (HDA).⁴⁴ The HDA scheme can be employed in both the traditional Casida approach and the polTDDFT approach, in the latter case by simply correcting the differences between KS orbital energies by the non-local exchange term. In short, HDA allows one to run B3LYP calculations with STO basis sets with a standard CPU usage using both Casida and polTDDFT algorithms. In order to facilitate the comparison with respect to the experiment, the Casida discrete spectra have been broadened with Lorentzian functions of 0.150 eV of Full Width at Half Maximum (FWHM). The polTDDFT results are intrinsically broadened by a Lorentzian function, and in order to have the same 0.150 eV FWHM, we employed in Eq. (5) an imaginary frequency ω_i of 0.075 eV.

Figure 1 shows the photoabsorption spectra of [Ag₂₄Au(DMBT)₁₈] anion, simulated via TDDFT using the Casida formalism⁵⁵ and the LB94⁵³ and B3LYP^51,52 xc-functionals in both the Kohn–Sham equations and the TDDFT kernel, together with the experimental spectrum taken at room temperature.⁴² In panels (a) and (b), discrete bars are shown together with the continuous spectrum obtained by broadening these peaks (with a FWHM of 0.15 eV), while in box (c) the continuous spectra simulated using both xc-functionals are compared with the experimental profile. Note that the theoretical spectra are simulated using a fully DFT-relaxed geometry of the MPC anion. The LB94 spectrum in panel (a) was extracted including the 500 lowest eigenvalues of the Casida matrix, thus spanning an excitation energy interval up to 3.52 eV (wavelengths >352 nm). In this case, it is not numerically safe to extract more than 500 eigenvalues of the Casida matrix because numerical instabilities start to appear. The B3LYP spectrum in panel (b) was obtained including only the 20 lowest eigenvalues of the Casida matrix, thus reaching an excitation energy of 2.46 eV (504 nm), before incurring into excessively long calculations. The reason for such a drastic reduction of the spectral energy interval accessible using B3LYP with respect to LB94 is connected with issues when using a hybrid xc-functional in our codes. When using a hybrid xc-functional, in fact, the efficiency of the ADF code is strongly hampered by the use of Slater-type orbital (STO) basis functions, which require a computationally heavy fitting procedure to calculate the non-local exchange in the response kernel. Fortunately, we have recently shown that this limitation can be overcome by employing the Hybrid Diagonal Approximation,⁴⁴ as we will see hereafter. In panel (c), the simulated broadened photoabsorption profiles are compared with experiment.⁴² Although the energy span of the B3LYP spectrum is quite limited, we are nevertheless able to assess the accuracy of TDDFT using two different xc-functionals on the lowest-energy peak, which in experiment falls at 619 nm (2.00 eV), while theory predicts it at 765 nm (1.62 eV) at the LB94 level and at 685 nm (1.81 eV) at the B3LYP level, respectively. It is clear that for this spectral feature, B3LYP performs much better than LB94 when compared with experiment, yielding an error of ∼0.2 eV, which is half the error given by LB94. Quite interestingly, the prominent experimental feature at 468 nm (2.65 eV) is very nicely reproduced by LB94, which predicts it at 459 nm (2.70 eV). The B3LYP simulations unfortunately do not reach this value of excitation energy, and therefore, it is not possible to compare B3LYP predictions for this peak. From this preliminary analysis, we conclude that the low energy features are reproduced reasonably well by the hybrid B3LYP xc-functional but not by the local density approximation nor even Coulomb-tail-corrected, the LB94 xc-functional. LB94, however, appears to perform well for the features at higher energy, as noted in previous studies (see Ref. 38 and references quoted therein). What, to the best of our knowledge, cannot be found in previous literature and is still missing is: (1) to extend this analysis to higher excitation energies; (2) to analyze the different nature of the transitions involved to understand why the low energy spectral feature is not well described at the non-hybrid level. Hereafter, we will try to give an answer to these open questions.

As a next step, we investigated the effect of the choice of the cluster geometry on the photo-absorption spectrum. Our previous analysis³⁸ showed that, in the case of atomically precise nanoclusters made of Au and aliphatic ligands, in order to obtain a good agreement between simulations and experiment in the calculation of photoabsorption spectra, it is important to employ the experimental geometry in the calculations (at least for the cluster core, i.e., the metal atoms and the sulfur atoms in the ‘staple’ fragments), rather than a fully DFT-relaxed geometry. Indeed, a worse agreement between theoretical and experimental spectra was obtained when fully optimized geometries were employed in the simulations in Ref. 38. Figure 2(a) shows the photo-absorption spectra, simulated using the LB94 xc-functional on the fully optimized geometry and on a partially optimized geometry, in which the cluster core (Ag, Au, and S atoms) is kept fixed and only the organic residues are optimized, together with the experimental spectrum taken at room temperature. To avoid the limitations connected with the diagonalization of large Casida matrices, whose difficulties hamper the calculation of the spectrum beyond 500 eigenvalues for LB94, and thus prevent access to features above 3.52 eV, in Fig. 2 and hereafter, we calculated all the spectra using the polTDDFT method, which does not have energy limitations and whose accuracy has been validated to be within ∼0.1 eV in the peak position with respect to the Casida approach.⁵⁷ In this connection, it must be noted that the accuracy of polTDDFT is now improved with respect to the original tests of the approach because it benefits from the use of the new auxiliary basis sets (here employed in Figs. 4, 9, and 10). To illustrate this point, we underline that polTDDFT using the new auxiliary basis set now provides an essentially quantitative agreement with respect to the Casida results even for very small gold clusters, such as, for example, Au₆ (as shown in Fig. S3 of the supplementary material), whereas the polTDDFT error in peak positions was up to 0.5 eV for the Au₂ system when using the original auxiliary basis set.⁵⁶ As apparent in Fig. 2, the agreement between the calculated and measured spectra improves when using the experimental geometry in simulations. We thus conclude that, also in the present case of an alloyed MPC and conjugated ligands analogously to what has been observed previously for pure Au MPCs and aliphatic ligands,³⁸ a fixed-core partially optimized geometry gives a better agreement with experiment, especially for the lowest-energy spectral feature, which is very sensitive to the choice of the geometry and is also the most difficult to describe by theory. In detail, using polTDDFT/LB94, the position of this low-energy peak is 1.60 eV for the fully optimized geometry, 1.73 eV for the fixed-core partially optimized and 2.00 eV at the experimental level. Note that for the prominent feature at higher energy (around 2.65 eV), we still find a very good agreement between theory and experiment, so this agreement is not sensitive to the choice of geometry, the differences being all below 0.1 eV. Regarding even higher excitation energies, we find that the smooth experimental structure around 3.75 eV has counterparts in the calculated profiles (in the range between 3.30 and 3.70 eV), whereas the sudden experimental increase with a steep slope above 4.2 eV is not properly reproduced by theory. In Fig. 2(b), we show the spectrum also in the wavelength scale, and in this format, the discrepancy between LB94 predictions and experiment for the low-energy (high-wavelength) feature is even more evident. We have also checked if the same conclusion, i.e., a better agreement between theory and experiment is obtained when using the experimental vs the DFT-optimized geometry, holds for TDDFT/B3LYP simulations as it does for TDDFT/LB94 simulations. The comparison of TDDFT/B3LYP results using the original and the new auxiliary basis set is reported in Fig. S1 of the supplementary material and fully supports the conclusion drawn from the LB94 analysis.

Having assessed that the experimentally constrained geometry gives a better agreement with experiment, in Fig. 3 we then investigate the effects of the choice of the xc-functional, while keeping frozen the experimental geometry of the cluster core and optimizing only the organic residues as described in Sec. III—panel (a) refers to the energy plot, and panel (b) refers to the wavelength plot. Besides the LB94 xc-functional, we considered the standard GGA PBE functional⁵⁰ as well as the hybrid B3LYP^51,52 within the HDA:⁴⁴ as noted at the beginning of this section, employing the HDA together with the polTDDFT scheme allows us to study the complete energy range of interest, overcoming the severe limitations of the standard Casida B3LYP scheme in the accessible excitation energy range. From Fig. 3, it is apparent that LB94 and PBE profiles are very similar, with a large error for the low energy feature with respect to experiment. On the contrary, the polTDDFT/B3LYP-HDA profile performs quite nicely with respect to experiment on the complete energy range. The low energy structure is calculated at 2.04 eV by B3LYP-HDA, in excellent agreement with experiment (2.00 eV). Moreover, not only the other features at higher energy remain in nice agreement with the experiment, but one can also see an increase in the oscillator strength intensity at 4.5 eV that parallels the sudden slope increase above 4.2 eV in experiment, although with a blue shift of about 0.3 eV to higher energy. The excellent result obtained with the B3LYP-HDA scheme is particularly evident when photo-absorption is plotted against wavelength, as shown in the lower panel of Fig. 3(b).

For completeness, in Fig. 4, we investigate the role of the different approximations to TDDFT/B3LYP, and specifically of B3LYP-HDA, within the ADF code, to identify precisely which scheme affords the best compromise between accuracy and efficiency. First, we considered the pure effect of using the HDA, i.e., still within a Casida scheme (not polTDDFT): this scheme predicts the lowest-energy spectral feature at 1.96 eV (632 nm) vs 1.81 eV (685 nm) at the full B3LYP Casida (without HDA, see Fig. 1). The HDA per se thus introduces an error of only 0.15 eV, in line with previous work.⁴⁴ Keeping now the HDA but employing the polTDDFT scheme instead of the Casida one, the low energy feature is calculated at 2.04 eV (608 nm), so the polTDDFT approximation introduces a discrepancy of only 0.08 eV with respect to the Casida HDA scheme. This result was obtained using the original fitting basis functions derived in Ref. 44. We have recently derived a new set of density fitting basis functions purposely optimized for polTDDFT/B3LYP-HDA calculations, a brief description of this set can be found in the supplementary material. Further discussion of these new fitting functions will be the subject of a future publication but details can also be found in the code documentation as they have been included in the 2021 version of the AMS [https://www.scm.com/product/ams/] suite of programs.^65,66 The polTDDFT/B3LYP-HDA profile obtained using this new set of density-fitting auxiliary functions (green curve in Fig. 4) is improved with respect to the old set and is in excellent agreement with the TDDFT/B3LYP-HDA Casida profile, giving a low-energy feature at 1.94 eV (639 nm) with an error of only 0.02 eV with respect to experiment together with a significant improvement in the prediction of the absolute intensity of this and other peaks.

Having so derived a very robust computational protocol, it remains as a scientifically important question to identify the nature of the most relevant spectral features in the photoabsorption spectrum. Moreover, it is also very important to understand why the low energy feature is the most difficult to predict via theory. A summary of the transition assignments, from the polTDDFT/B3LYP/HDA calculations with the new density fitting basis set, are gathered in Fig. 5. The details of the transitions are provided in Table I. The transition at 1.94 eV is ascribed to a mixed character HOMO → LUMO + 1, HOMO − 1 → LUMO, HOMO − 2 → LUMO, and HOMO − 1 → LUMO + 1. The B3LYP orbitals involved in this and the following transitions are depicted in Fig. 6 (occupied ones) and Fig. 7 (virtual ones), respectively. They can be characterized as follows: the occupied orbitals are strongly mixed with high Ag and S contributions and can be classified as Ag–S bonds at the staple units with the involvement of Ag atoms belonging to the cluster. The involved unoccupied levels are localized in the same spatial region of the occupied ones. The transition at 2.20 eV is weak and is not detected in the experiment. The next transition at 2.71 eV can be ascribed to HOMO − 1 → LUMO + 3 and HOMO → LUMO + 4. Interestingly, the involved occupied orbitals are the same as in the low-energy peak, whereas the virtual orbitals are completely different. Note that the LUMO and LUMO + 1 are almost degenerate, whereas the LUMOs at higher energy are well separated. From the plot of the virtual orbitals, it is apparent that the LUMO + 3 and LUMO + 4 still have strong Ag and S contribution, but they are shifted toward the cluster core with respect to the staple units. The next transition at 3.38 eV is instead ascribed to π → π* ligand transitions, with the virtual π* strongly mixed with S–C bonds and Ag atoms, so that for such transition some ‘charge transfer’ character from ligands to metal is found. As in the first three absorption bands the same occupied orbitals are involved in the transitions, we can conceivably ascribe the difficulty in describing the low energy bands via TDDFT/LB94 to an error in the energy positions of the virtual orbitals. LB94 is indeed known to excessively stabilize the lowest virtual orbitals (LUMO and LUMO + 1) and to place them at too low energies, whereas B3LYP, containing a Hartree–Fock non-local exchange component, raises the energy of virtual orbitals thus yielding more accurate results. We can further conjecture that such effects must compensate each other for the orbitals beyond LUMO + 2 since the higher experimental features are properly described by both LB94 and B3LYP. However, to identify, with precision, the reasons for the better performances of B3LYP in comparison with standard or tail-corrected GGA approaches (such as PBE and LB94), we show in Table II the outcome of an analysis of the low energy transition predicted by the LB94 level at 1.74 eV. Interestingly, several similar excited configurations are involved in this transition, although with different weights with respect to B3LYP. More surprisingly, a comparison between Figs. 6 and 7 with Fig. 8 (where we show the orbitals calculated at the LB94 level) shows that the virtual orbitals are almost identical from the two methods, whereas the occupied ones (HOMO, HOMO − 1, and HOMO − 2) are mainly S 3p in LB94 and Ag 4s in B3LYP. Our conclusion is that the B3LYP functional is more accurate in describing the electronic structure of these clusters, with a more balanced contribution from Ag and S, while LB94 tends to overestimate the role of S 3p in the HOMO and HOMO − 1 orbitals by putting them at too low energies. This is the origin of the too low energy of the LB94 excitations at the absorption threshold, which are instead correctly placed by B3LYP in the proper (higher) energy region, overlapping with other excitations involving the same sets of occupied and virtual orbitals. Clearly, this difference has consequences also on the coupling between the orbitals, which takes place at the TDDFT level, and indeed the composition of the transition in terms of one-electron excited configurations is also somewhat changed in going from LB94 to B3LYP (Hartree–Fock exchange terms in the diagonal matrix elements of the TDDFT kernel change linear response both affecting directly excitation energies and indirectly by changing the composition of the excited states⁴⁴). Therefore, we are led to ascribe the better performances of B3LYP to its more accurate capacity to describe the electronic structure of the clusters, in particular the balance between metal and ligand contributions.

To complete our validation study, since the TDDFT-HDA/B3LYP scheme has proven the best choice for the present mixed and conjugated MPC compounds, we decided to check its performances also for the [Au₂₅(SR)₁₈]⁻ system, a non-conjugated and pure Au case that was investigated in our previous work.³⁸ In Fig. 9, the experimental data obtained for the [Au₂₅(SCH₂CH₂C₆H₅)₁₈]⁻ TOB⁺ system in the present work are compared with the simulated spectrum of the simpler [Au₂₅(SCH₃)₁₈]⁻ anion, using an experimental geometry for the AuSC core⁶⁸ and optimizing only the positions of the hydrogen atoms. In fact, in our previous work,³⁸ we have shown that for this aliphatic ligand (the phenyl ring is not directly bonded to the sulfur atom, but is tethered through an ethylene bridge), it is reasonably accurate to model the ligand by keeping the experimental geometry but simplifying it with the smaller SCH₃ fragment. This does not drastically deteriorate the appearance of the spectrum at low energy. From Fig. 9, it is apparent that the polTDDFT B3LYP HDA results with the new basis follow very accurately the Casida profile (more accurately than the old basis). The comparison with the experiment is also fairly nice, with the first two experimental features detected at 1.80 and 2.75 eV well reproduced by theory at 1.82 and 2.62 eV, respectively. The next experimental feature at 3.1 eV can be ascribed to the calculated structure in the range between 2.84 and 2.94 eV, obtaining quantitative match at low energy and a discrepancy about 0.2 eV at this higher energy.

As a final remark, it is important to note that calculations give spectral features that are more resolved with respect to the experiment. This is due to the fact that simulations correspond to frozen geometries or formally a temperature of 0 K, whereas experiments are conducted at room temperature and are therefore affected by thermal broadening of spectral features. For a more accurate comparison and validation of theory vs experiment, it is important to perform photo-absorption experiments at low temperatures.

We have therefore performed a series of experiments aimed at measuring low-temperature photoabsorption spectra of the Ag₂₄Au(DMBT)₁₈ single negative anion (experimental details are provided in the supplementary material). We collected spectra down to 77 K, and a series of representative spectra obtained at different temperatures is shown in Fig. S2. The comparison between the absorption spectrum experimentally measured at 77 K and simulated via polTDDFT/B3LYP using the new auxiliary basis set at the core-frozen geometry is reported in Fig. 10. This comparison shows an excellent accuracy, which is extremely encouraging because the much better resolved spectral features observed in the experimental spectrum compare even better with the features of the simulated spectrum. As the only discrepancy, a composite feature of the low-energy peak may be noted, exhibiting a small feature at around 740 nm, which is not present in the simulated spectrum. This discrepancy may be due to three main reasons: first, geometry effects connected with the difference between the solid-state and the solution environment. Indeed, the experimental spectra in Fig. S2 clearly show how sensitive is the lowest-energy peak with respect to vibronic coupling. Second, our predictions are based on a Lorentzian broadening of the theoretical spectra with a constant energy value throughout the excitation interval: especially at low energies, this may hide fine details and features when the energy splitting is below the broadening. Finally, Spin–Orbit Coupling (SOC) effects could also cause split peaks especially in the low-energy region.⁶⁸ At present, we cannot include SOC effects into our modeling, but we are currently working to solve this issue.

Atomically precise nanoclusters, or monolayer-protected clusters (MPCs), represent an ideal playground to assess and validate the performance of computational methods via a stringent comparison with chemical–physical experimental data, enabling singling out accurate simulation protocols that can then be applied with predictive confidence in the rational design of the optical properties of new metal-cluster-based materials.^1–39 Here, we focus on the prediction of the optical response (the photo-absorption spectrum) of these systems, we consider a prototypical compound, the Ag₂₄Au(DMBT)₁₈ single negative anion (DMBT = 2,4-dimethylbenzenethiolate ligand) as an alloyed and ligand-conjugated MPC, and we conduct a validation study to investigate which DFT/TDDFT simulation tools provide the best agreement between predicted and measured spectra.

For this purpose, we use efficient simulation algorithms, i.e., the complex polarizability polTDDFT approach and the Hybrid-Diagonal Approximation (HDA), that allow us to employ a variety of exchange-correlation (xc-) functionals with the needed accuracy and an affordable computational cost. We are so able to demonstrate that the spectrum of this prototypical compound, especially but not exclusively in the absorption threshold (low-energy) region, is sensitive to the choice of both the xc-functional and the MPC geometry. Technically, we confirm that the polTDDFT and HDA approaches achieve an accuracy within ≈0.1–0.2 eV in the position of the main photo-absorption peaks that are accessible to an exact TDDFT prediction.

Moreover and importantly, the comparison of simulated and experimentally determined spectra recorded from room temperature down to 77 K allows us to validate that a hybrid xc-functional employed in both the Kohn–Sham equations and the diagonal TDDFT kernel at the crystallographically determined experimental geometry (as proposed in Refs. 37 and 38) is able to provide a consistent agreement (∼0.1 eV accuracy) between simulated and measured spectra across the entire optical region. This level of accuracy is also obtained for a Au₂₅(CH₃)₁₈ single negative anion, also included in this work as an example of a pure Au and non-ligand-conjugated MPC.

Finally, we conduct a single-particle decomposition analysis of the results and are so able to trace the physical reasons for the failure of non-hybrid xc-functionals in predicting the optical response of these systems. These are found to be due to the inability of non-hybrid xc-functionals to correctly describe the balance between metal and sulfur contributions to the electronic structure of these systems, and in particular to their tendency to put S 3p orbital contributions to bands around the HOMO–LUMO gap at too low energies, thus underestimating the position of the corresponding excitations that are instead correctly blue-shifted and merged with higher-energy peaks when hybrid xc-functionals are employed.

See the supplementary material for tests using the B3LYP xc-functional with fully optimized geometry, Cartesian coordinates of the structures used in the TDDFT simulations, numerical raw data for Figs. 1–4 and 9, and brief description of the new fitting set.

Computational support from the CINECA supercomputing center within the ISCRA program is gratefully acknowledged. The authors are grateful to the Stiftung Beneficentia for a generous grant employed to set up a computational server. Support from the University of Trieste within the FRA program is gratefully acknowledged.

The data that support the findings of this study are available from the corresponding authors upon reasonable request.