Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX-XXX
pubs.acs.org/JCTC
Optical Gaps in Pristine and Heavily Doped Silicon Nanocrystals: DFT versus Quantum Monte Carlo Benchmarks R. Derian,† K. Tokár,† B. Somogyi,‡ Á . Gali,‡ and I. Štich*,†,§ †
Center for Computational Materials Science, Institute of Physics, Slovak Academy of Sciences, 84511 Bratislava, Slovakia Wigner Research Centre for Physics, Institute for Solid State Physics and Optics, Hungarian Academy of Sciences, Budapest, Hungary § Department of Natural Sciences, University of Ss. Cyril and Methodius, 91701 Trnava, Slovakia ‡
S Supporting Information *
ABSTRACT: We present a time-dependent density functional theory (TDDFT) study of the optical gaps of lightemitting nanomaterials, namely, pristine and heavily B- and Pcodoped silicon crystalline nanoparticles. Twenty DFT exchange−correlation functionals sampled from the best currently available inventory such as hybrids and rangeseparated hybrids are benchmarked against ultra-accurate quantum Monte Carlo results on small model Si nanocrystals. Overall, the range-separated hybrids are found to perform best. The quality of the DFT gaps is correlated with the deviation from Koopmans’ theorem as a possible quality guide. In addition to providing a generic test of the ability of TDDFT to describe optical properties of silicon crystalline nanoparticles, the results also open up a route to benchmark-quality DFT studies of nanoparticle sizes approaching those studied experimentally.
1. INTRODUCTION In the last 20 years, porous Si and Si nanocrystals (Si-NCs) have played an important role as environmentally friendly lightemitting nanomaterials.1 Si-NCs may exist in a variety of forms, with colloidal Si-NCs being perhaps the most widely studied form, with a diversity of applications including use as a precursor for Si-NC films for optoelectronic applications. By modification of the growth temperature, the diameter of the SiNC particles (Si-NCPs) can be varied over a wide range of ∼2−15 nm.2 For applications, Si-NCPs are heavily B- and P-codoped. The doping, among other effects, helps to protect them in solution and modifies their photoluminiscence (PL) properties.1 For instance, the codoped colloidal Si-NCs exhibit size-controllable PL over a very wide energy range of 0.85−1.85 eV,1 i.e., including PL well below the bulk Si band gap of 1.12 eV. Since quantum confinement can only increase the band gap,3 the below-bulk-band-gap PL provides evidence that the PL arises from transitions between the donor (P) and acceptor (B) states in the band gap of Si-NCs.1 The donor/acceptor localization of the charge of the HOMO/LUMO states means that these transitions are of the charge-transfer (CT) type (see below for more details). Hence, in addition to the nanoparticle size, the PL also depends on the relative B:P:Si concentration ratio4 and dopant distribution5 in the NCPs. Furthermore, the PL properties may also be affected by NC-to-NC interactions causing exciton migration between NCs by Förster energy transfer.6 © XXXX American Chemical Society
Because of all these complexities, atomic-scale understanding of the properties of Si-NCPs by experimental methods1 would greatly benefit from theoretical modeling. The nanometer size of the NCPs represents a hurdle for such modeling. The smallest experimental NCs have diameters of ∼2 nm, and the usual size is larger than 4 nm;1 hence, ∼4000−20000 electrons have to be handled. These sizes can be computationally approached using time-dependent density functional theory (TDDFT).7 However, description of band gaps using DFT may be heavily biased. For example, the SC40 data set8 suggests mean absolute gap errors of commonly used hybrid DFT exchange−correlation (xc) functionals between 0.37 (HSE9) and 0.59 eV (PBE010), indicating on one hand a strong dependence on the xc functional used but on the other hand the possibility of proper tuning of the functional via results of more accurate modeling, such as GW11 or quantum Monte Carlo (QMC).12 While more accurate than DFT, the GW gaps in customary approximations may also be strongly dependent on the underlying DFT xc functional and on other approximations.11 The QMC method has a documented ability to capture both static and dynamic correlations on both ground and excited states, which are often important for accurate gap values.12−15 This type of DFT comparison/tuning versus benchmark-quality theory has, compared with theory−experiment comparison/tuning, the advantage that the same approximations, such as the vertical transition approximation, Received: August 4, 2017 Published: October 30, 2017 A
DOI: 10.1021/acs.jctc.7b00823 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
2. COMPUTATIONAL DETAILS With the quest to benchmark the performance of DFT xc functionals to describe excited states, the optical gaps in particular, in pristine and codoped Si-CNPs against robust QMC results, we selected perhaps the smallest reasonable structures for study: the Si38H42/Si36BPH42 cluster models (Figure 1). The structures optimized at the PBE level24,25 were used in all of the electronic structure calculations.
absence of solvent, vibronic effects, etc., can be made in both theoretical descriptions. The main limitation of the QMC method is its computational cost, which limits the size of the model system to just a few hundred electrons. While TDDFT has in the past decade reached its maturity and several test sets have been developed,16,17 none of them appears to cover the more “solid-state”systems of interest here. Studies of excited states on these sets suggest that the hybrid functionals improve the description of excited states, including the electronic gaps, compared with the pure functionals and that the hybrid functionals outperform the range-separated hybrids.16 On the other hand, the range-separated hybrids may help to describe the CT excited states.16,18 Most of these previous studies dealt with singlet−singlet transitions and more recently also with singlet−triplet transitions.19 Experimentally, both singlet−singlet and singlet−triplet transitions can be probed,14,20 and therefore, here we study both transitions. The functional dependence is significantly worse with TDDFT than with DFT as a result of the adiabatic approximation used in TDDFT.21 Specifically, with hybrids there is a strong dependence on the amount of exact exchange (EE) mixed in:22 small-gap systems with efficient screening may require less exact exchange, and on the contrary, large-gap systems with inefficient screening may require more than the canonical ∼20% EE. Such a situation will occur also in the codoped SiNCs, where the gaps, depending on the size of the nanoparticle, vary by more than 1 eV, reflecting the change in the dielectric screening.1 In addition, in the context of higher-lying states, such as Rydberg states, the basis set dependence on diffuse functions has been emphasized.23 The paper has the threefold purpose of (1) finding the best xc functional(s) for a truly large-scale treatment of Si-NCPs via benchmarking against ultra-accurate QMC results on model systems, (2) correlating these findings with alternative cheaper indicators of the quality of results, and (3) providing another specialized test set for DFT functionals specifically targeting optical gaps. To this effect, we explicitly calculate the optical gaps of small model pristine and B- and P-codoped Si-NCPs with a number of DFT xc functionals and benchmark the results against the QMC results. The selected xc functionals bracket the most important commonly used families, such as local (L), generalized gradient approximation (GGA), metaGGA (M), hybrid (H), and range-separated hybrid (RSH). In view of the very large number of available xc functionals, of particular interest will be comparison of different members of the same family, such as hybrids or range-separated hybrids, to see whether they yield similar quality of results. In a situation where many dozens of xc functionals with vastly different properties are available, benchmarking the DFT xc functionals against more robust techniques such as correlated quantum chemistry, GW, or QMC should be routinely done. However, given their prohibitive cost, simpler quality descriptors based purely in the DFT domain would be desirable. Of particular relevance for optical gaps is the relation between the ionization threshold (−ϵDFT HOMO) and ionization potential (IP), where in DFT −ϵDFT HOMO often is much lower than IP, causing gaps to close at too-low values. This relation as a quality indicator of the optical gap description will be scrutinized. More sophisticated procedures that are similar in spirit have recently been used to tune xc functionals to describe CT transitions.18 Furthermore, extrapolation of our results obtained on small model silicon crystalline nanoparticles to larger systems, including those encountered experimentally, will be discussed.
Figure 1. Structural model of a H-terminated B- and P-codoped model Si-CNP, Si34BPH42. The B and P atoms are shown in green and gray, respectively.
To set the stage, we first determined the nature of the HOMO → LUMO transition. Figure 2 indicates that the transition is of the CT type, as can clearly be seen from the charge localization primarily around the acceptor/donor sites for the HOMO/LUMO states, respectively. The optical gaps were calculated using the TDDFT technique.7 A few tests were also performed making use of the Tamm−Dancoff approximation (TDA).26 TDA was found to improve the singlet−triplet gaps in molecular systems.26 Our model nanoparticles represent more a “solid-state” context, and therefore, we wanted to test whether a similar conclusion also holds for our systems. The lowest 20 singlet/triplet excited states were optimized for the pristine and B- and P-codoped nanoparticles. For benchmarking purposes, we selected a wide range of DFT xc functionals from different major families. For clarity, in addition to their names, the xc functionals are labeled by the color-coded type and the xc functional number. The local functional (L, olive green) is SVWN5 (#1);27 the GGA functionals (GGA, orange) are PBE (#2),24,25 PW91 (#3),28 and BP86 (#4);29,30 the meta fuctionals (M, sky blue) are TPSS (#5),31 VSXC (#6),32 and M06-L (#7);33 the hybrid functionals (H, lime green) are B3LYP (#8),34,35 X3LYP (#9),36 B98 (#10),37,38 APFD (#11),39 PBE0 (#12),10,40 M06 (#13),41 BMK (#14),42 and M06-2X (#15);41 and the range-separated hybrids (RSH, pink) are HSE06 (#16),9 ωB97X-D (#17),43 M11 (#18),44 CAM-B3LYP (#19),45 and LC-ωPBE (#20).46 The functionals are ordered by increasing fraction of EE included (H) or increasing ω parameter defining the range separation (RSH). By construction, these xc functionals have different treatments of electron correlation. However, the largest single difference lies in the treatment of electron exchange. Consequently, they have different degrees of selfinteraction error,47,48 signaled also by collapse of Koopmans’ B
DOI: 10.1021/acs.jctc.7b00823 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
Our TDDFT and QMC gap estimates both use the vertical transition approximation, i.e., they neglect all adiabatic, vibronic, and zero-point vibrational energy effects. Such approximations tend to increase the gap value compared with experiments, in some cases quite significantly.16,17 However, since we are consistently comparing gaps calculated within the vertical transition approximation, neglecting the vibronic effects does not affect our conclusions. For CI modeling we used the GAMESS suite of codes,50,51 while for all of the VMC and DMC calculations the QWalk code53 was used. In all of the calculations, the atomic cores were replaced by effective core potentials (ECPs)54 used in combination with VTZ basis sets. The truncated CI expansion led to a single CSF for the ground-state wave function, S0, whereas the excited-state expansion led to four CSFs for the singlet, S1, and two CSFs for the triplet, T1. Dynamical correlations were explicitly built into the trial wave functions via the isotropic Schmidt−Moskowitz Jastrow factor,52,55 including electron−electron, electron−nucleus, and electron−electron− nucleus correlations. The parameters of the trial wave functions,52 including the weights of the CSFs, were optimized by minimizing linear combination of energy and variance.56 The T-moves scheme57 was used to keep the calculations with ECPs variational.
3. RESULTS AND DISCUSSION The performance of TDDFT with the different xc functionals in describing the singlet and triplet optical gaps is shown in Figures 3 and 4. More details, including the TDA results, densities of excited states and optical spectra, are provided in the SI. As discussed above, doping results in reduction of the electronic gap. This feature is correctly reproduced by all descriptions. In most cases the reduction is by ∼0.8 eV, compared to the DMC reduction by ∼0.9 eV. With very few exceptions, the trends in the performance of the xc functionals are qualitatively similar for both pristine and codoped particles as well as for singlet and triplet gaps. On average, the errors are larger for the singlet gaps than for the triplet gaps and for the pristine CNP than for the codoped one (see Figures 3 and 4). At variance with molecular test sets,17 which uniformly favor the hybrids or double hybrids, we find the RSH functionals to perform best. Indeed, the errors are largely reduced in the RSH group compared with the H group, with the former slightly overestimating the gap and the latter hugely underestimating it. The exception among RSHs is HSE06 (#16), which underestimates the gaps similarly to the hybrids. This distinction is due to its construction, which puts the exact exchange into the short-range part rather than the long-range part, which was designed to be similar to PBE0 while reducing computational cost.9 These findings agree with previous conclusions for CT transitions.16,18 Overall CAM-B3LYP (RSH, #19) is found to be the best performing xc functional, followed by ωB97X-D (RSH, #17), M06-2X (H, #15), and BMK (H, #14) for singlet gaps only. This means that the best H functionals for our test set are those with a very large EE portion (54% in M06-2X and 42% in BMK). This is in line with the previous conclusions that CT transitions may require xc functionals with a large portion of EE.16 However, unexpectedly and at compete variance with previous results that favor the BMK xc funtional (H, #14) for singlet−triplet transitions,19 we find the BMK functional to be less satisfactory for the singlet−triplet transition, especially for the codoped nanoparticle. While all of the GGA and meta functionals have very similar performance, the RSHs have more
Figure 2. CT type of the HOMO → LUMO transition illustrated by CI natural orbitals: (A) the HOMO is localized primarily on the acceptor (B site); (B) the LUMO is localized primarily on the donor (P site).
theorem,49 which we carefully monitored. In addition, we also add results of time-dependent Hartree−Fock (TDHF) and configuration interaction with single and double excitations (CISD) calculations. Hence, in total we calculated 80 excited states (pristine, codoped, singlet, and triplet (20 states each), all calculated with 20 different xc functionals). In the main text, we benchmark only the electronic gaps; the densities of excited states and the corresponding optical spectra are shown in the Supporting Information (SI). For CI, DFT, and TDDFT/ TDHF modeling, we used the GAMESS suite of codes.50,51 The fixed-node QMC12,52 electronic gap benchmarks were performed in three steps: (1) in order to capture also the possible static correlations and multireference character, the trial wave functions for the ground and excited states were constructed from CISD expansions, keeping all of the configuration state functions (CSFs) with weights larger than 0.05; (2) the trial wave function constructed from a truncated CISD expansion was optimized using variational Monte Carlo (VMC) techniques; and finally, (3) optical gaps were computed from the diffusion Monte Carlo (DMC) energies of the ground and excited states as the first singlet/triplet vertical excitation energies: s/t Δopt ≈ Evs/t = E1s/t − E0s
C
DOI: 10.1021/acs.jctc.7b00823 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation
Figure 3. Errors with respect to the DMC value for the TDDFT singlet optical gaps of the pristine (upper panel) and codoped (lower panel) model Si-CNPs calculated with different xc functionals. For completeness, TDHF and CISD results are also shown. Results for the different xc functionals are color-coded (see the inset) and numbercoded underneath the panel. The functionals are ordered by fraction of exact exchange included (H) or by the ω parameter defining the range separation (RSH).
Figure 4. Errors with respect to the DMC value for the TDDFT triplet optical gaps of the pristine (upper panel) and codoped (lower panel) model Si-CNPs calculated with different xc functionals. For completeness, TDHF and CISD results are also shown. Results for the different xc functionals are color-coded (see the inset) and number-coded underneath the panel. The functionals are ordered by fraction of exact exchange included (H) or by the ω parameter defining the range separation (RSH).
mixed performance. In all cases, CAM-B3LYP (RSH, #19) and ωB97X-D (RSH, #17) have consistently better performance than the other three RSH functionals studied, HSE06 (RSH, #16), M11 (RSH, #18), and LC-ωPBE (RSH, #20). The observed differences can be traced back to the amount of EE incorporated via the value of ω used in construction of the RSH functional, keeping in mind that the functionals may start with radically different amounts of EE, such as M11 with 40% and the HSE06/PBE0 similarity by construction. For example, LCωPBE overestimates the excitation energy and ωB97X-D lowers the trend as a result of a smaller ω value and less EE included, confirming again that it is primarily the amount of EE incorporated that dictates the quality of the results. Perhaps surprisingly, the performance of TDHF very similar to or better than that of GGA and meta-GGA TDDFT (also see the discussion of Koopmans’ theorem below). Hence, while our general conclusion of the preference of RSHs for our CT-type of transition is similar to those found previously, other findings, such as the clearly best performance of CAM-B3LYP and ωB97X-D and the poorer performance of the BMK for triplets, makes a priori selection of an xc functional for gap calculation based purely on its construction type impossible. While the L, GGA, and M functionals are expected to show poor
performance for CT-type gap calculations, our benchmarks indicate that the subtle details in the treatment of the exchange, which is the leading term in the gap formation, are crucial for performance of the H and RSH xc functionals. Perhaps surprisingly, TDA26 improves the gap values only negligibly; the largest gap correction was
