Dependence of Nonadiabatic Couplings with KohnâSham Orbitals on

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Dependence of Nonadiabatic Couplings with Kohn-Sham Orbitals on the Choice of Density Functional: Pure vs. Hybrid Yuhan Lin, and Alexey V. Akimov J. Phys. Chem. A, Just Accepted Manuscript • Publication Date (Web): 31 Oct 2016 Downloaded from http://pubs.acs.org on October 31, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Dependence of Nonadiabatic Couplings with KohnSham Orbitals on the Choice of Density Functional: Pure vs. Hybrid Yuhan Lin and Alexey V. Akimov* Department of Chemistry, University at Buffalo, The State University of New York, Buffalo, NY 14260-3000 *Corresponding author. Email: [email protected]

Abstract

Nonadiabatic molecular dynamics (NA-MD) is an extremely useful approach to model electron transfer dynamics in molecular and solid state systems. The performance of NA-MD simulations depends critically on the accuracy of the underlying electronic structure properties, such as state energies and nonadiabatic couplings (NAC). Practical NA-MD modeling relies extensively on computationally efficient pure density functionals, despite the known intrinsic problems of the latter. A reliable solution to the problems presented by the use of pure density functionals, is the use of hybrid functionals, however hybrid functionals are significantly more expensive. In this work, we investigate how the choice of the density functional can affect the NAC magnitudes for small silicon clusters and silicon hydrides. We find that pure functionals can significantly 1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 36

overestimate NAC magnitudes in comparison to those obtained with hybrid functionals. The ratio of the two magnitudes increases with the system size and differs by an order of magnitude even for small Si7 and Si26 clusters. We present a detailed analysis of the correlations of the NACs computed with different methods, discuss the fundamental grounds for the observed differences, and propose a simple scaling technique for correcting NACs that can be used within the context of NA-MD simulations.

1. Introduction Sustainable energy production is among the main challenges of the 21st century due to the increasing energy demand, limited fossil fuel resources, and the amount of pollutants released when the latter are utilized. The conversion of solar energy into electricity by photovoltaic (PV) materials is an attractive route to address the global energy demand. In response to this idea, numerous experimental1–7 and theoretical1,8–14studies emerged as a result. A comprehensive understanding of the closely intertwined charge transfer (CT) and energy transfer (ET) processes in solar harvesting materials often requires explicit atomistic time-domain simulations. Nonadiabatic molecular dynamics (NA-MD)15–19 has been the method of choice for such simulations, describing the coupled evolution of electronic and nuclear degrees of freedom. The technique has found use in various applications such as the studying of excitation energy relaxation, charge separation in conjugated oligomers20–23 at interfaces,24–27 multiple exciton generation28–30 and Auger relaxation31–34 in quantum dots, singlet fission35–37 and charge transport38–40 in organic photovoltaic materials. Most NA-MD methods seek a detailed atomistic description, and depend heavily on the efficiency and accuracy of the underlying electronic structure calculations. The methodology employed varies broadly, depending on the particular application. High-level correlated ab initio 2 ACS Paragon Plus Environment

Page 3 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


methods for NA-MD applications have been extensively used by several groups, often in the context of photoinduced isomerization of organic molecules.41–53 Although the correlated multiconfiguration and multi-reference methods constitute a solid basis for the studies of excited state dynamics, their use is limited to relatively small systems, due to the rapidly scaling cost of computation. As an alternative to the expensive wavefunction-based methods, several groups employ simplified approaches based on the density functional theory (DFT) formulation. There are essentially two main branches on this front: a single-determinant Kohn-Sham (SDKS) approach developed in the Prezhdo group54–57 and the restricted open-shell Kohn-Sham (ROKS) approach initially invented by Frank and co-workers58 and later extended by others.19,59–61 An SDKS-type methodology was utilized by Shenvi and Tully62 in the context of simulating the dynamics of electrons in metals. For more extensive discussion of the DFT-based excited state dynamics methods, we refer the reader to an excellent recent review by Barbatti and Crespo-Otero.63 Although the SDKS and ROKS approaches are rather simplified, they have been particularly useful to model excited states dynamics in extended systems, provided the electronic excited states can be reasonably approximated by one (for SDSK) or two (for ROKS) Slater determinants .11,12,64–70 Apart from the computational efficiency, the SDKS and ROKS methods have the advantage of providing a well-defined and simple way of computing nonadiabatic couplings (NAC), which is needed for quantum dynamics simulations. In the simplest case, a numerical scheme of HammesSchiffer and Tully71 can be applied to compute scalar NACs (time-derivatives). It should be noted that NACs utilized in many DFT-based NA-MD simulations are computed in a simplified way. Namely, the electronic states are approximated as Slater determinants built of KS orbitals,

3 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

[

Page 4 of 36

]

∂ ΦJ , Φ I = Aˆ ϕi1 (1)ϕi2 (2 )...ϕi N ( N ) . The NACs computed with such wavefunctions, d IJ = Φ I ∂t can be reduced to NACs computed with KS orbitals, d IJ = ϕi

∂ ϕ j .54,56,72 KS orbitals appear ∂t

in DFT as auxiliary states of non-interacting electrons, so the physical meaning of these orbitals is undetermined. In many cases, these states resemble rigorously-defined wavefunction-based states (e.g. from the Hartree-Fock theory). The second notion is that the correlated wavefunction theories, like the configuration interaction (CI) or coupled cluster methods, require more than a single determinant to accurately describe many-body electronic wavefunctions. Yet, it is not uncommon that only a single configuration is dominant (has the highest weight in the CI expansion). Although rather simplistic, the above two approximations prove to be quite useful for practical studies, allowing many researchers to obtain new insights into nonadiabatic dynamics in many materials.11,12,64–68,73 Despite the common utilization of the NACs defined with KS orbitals as shown above, their potential dependence on the choice of density functional remains relatively unexplored. In the present work, we address one of the possible sides of this problem: whether the use of hybrid functionals significantly affects the magnitudes of the scalar NACs defined above and the general trends that exist among pure and hybrid functionals. A more extensive benchmarking using many-body wavefunctions can be envisioned, but that goes beyond the scope of the present work. Apart from possible strong correlation issues present in the SDKS and ROKS formulations due to their intrinsic approximations, there exists a question of dynamic correlation and its influence on the excited states dynamics. One of the main advantages of DFT over the existing postHartree-Fock methods is that the former provides an efficient way to account for electronic correlation via appropriately designed exchange-correlation functionals. The efficiency achieved


Page 5 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with the local density approximation (LDA) or generalized gradient approximation (GGA) functionals comes at the price of reduced accuracy, especially in regard to the prediction of excited states properties. The notorious problems of DFT with pure (LDA or GGA) functionals include the underestimation of the highest occupied molecular orbital (HOMO) - lowest unoccupied molecular orbital (LUMO) energy gaps74–78 and the related underestimation of the charge-transfer excited state energies.79–82 Significant improvements can be achieved by using hybrid functionals, which substitute a fixed (global hybrid) or geometry-dependent (rangeseparated hybrid) fraction of the DFT exchange with the exact (Hartree-Fock) exchange.76,83–85 For example, recent studies have demonstrated that hybrid functionals improve the excitation energies of CT states,86 are advantageous in computing sum-frequency generation87 and optical78 spectra, and help improve the size extensivity of DFT computations.88 Despite the advantageous properties of hybrid functionals, the computational costs they incur preclude their use in NA-MD calculations. In the most straightforward implementation, the costs of pure DFT calculations are determined by the linear algebra operations, leading to a cubic

( )

scaling with the number of orbitals, O N 3 . On the contrary, the costs for the Hartree-Fock

( )

exchange integrals scale quartically, O N 4 . Certain accelerations are possible, especially if the integrals are evaluated in the localized basis. For finite and moderately sized systems, only the short-range contribution of the hybrid functional will be dominant, while not adding significant computational overhead. As the size of the system increases, the long-range terms will come into play, adding the expected quartic dependence of computational costs. For the periodic systems and the delocalized plane-wave basis sets, the difference in costs of pure vs. hybrid functionals may quickly exceed an order of magnitude, since the range of exchange may well exceed the size of the unit cell. Thus, although less expensive variations are possible when hybrids are utilized 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 36

instead of pure functionals, the additional slow-down of the former in MD or NA-MD contexts makes these functionals less popular. Most calculations performed within the SDKS or ROKS formulations rely on pure DFT functionals. The associated problems with the underestimated energy gap can be resolved by introducing empirical corrections to electronic state energies.36,57 Within the trajectory surface hopping NA-MD formulations, such corrections affect only the phase evolution of the electronic states, contributing only to the pure dephasing. The electronic transition rates, modulated by the magnitude of NAC, remain practically unaffected. From the practical considerations, NACs are computed only at the GGA level, especially for extended systems. Computations of NACs using hybrid functionals require a significantly larger amount of time when compared to similar calculations performed with pure density functionals, especially when considering that thousands of MD steps are used routinely. At the same time, the magnitude of NACs computed with hybrid functionals may be different from those computed with pure functionals, hence affecting the rates of nonadiabatic transitions. The question of how large such difference could be remains largely unanswered. In the present work, we investigate the dependence of NACs computed using KS orbitals on the choice of exchange-correlation functional. We focus on the differences between NAC computations using pure and hybrid functionals. As a test set, we consider a number of silicon hydride molecules and nanoscale silicon clusters (quantum dots, QD). Silicon QDs constitute a promising class of materials for solar energy harvesting and light-emitting diode applications and hence are actively studied.89–94 Most importantly for this study, the size of the systems considered allows us to compute NACs thermally-averaged over sufficiently long molecular dynamics trajectories using hybrid functionals. Our study estimates the scaling coefficients


Page 7 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


needed to convert NACs computed at different levels of theory and establishes certain correlations. We also investigate how the scaling and correlations are affected by the chemical nature and size of the atomistic systems. 2. Methodology 2.1. Electronic Structure and Molecular Dynamics. In order to test the sensitivity of NACs to the choice of density functional, we consider a number of small silicon hydride molecules and nanoscale silicon clusters (Figure 1). Specifically, the following five systems are studied: SiH4, Si2H6, Si5H12, Si7, and Si26. The small size of the systems is essential, because our calculations involve hybrid functionals and long trajectories needed to obtain statistically-meaningful values. Moreover, our calculations are performed using a periodic code with a planewave basis. As mentioned in the introduction, the hybrid functional calculations may become expensive very quickly under such conditions. Hybrid functional calculations on larger systems are possible, but they are limited to either single point calculations or geometry optimizations, with only few nuclear steps, unless linear-scaling approaches are utilized. Despite a relatively small size of the systems, the electronic structure calculations with hybrid functionals are repeated in NAC calculations hundreds of times. Hence, the work represents a major computational effort needed to benchmark the performance of density functionals in their relationship to NA-MD simulations. The numerical results obtained for the present Si-based compounds are specific to these systems. It is reasonable to assume that the functional-dependent trends would persist for other systems as well, especially systems of a similar type. The conclusions of the present work are not expected to hold for significantly different chemistries, in which additional factors (e.g. relativistic effects, etc.) may play a dominant roles.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a)

Page 8 of 36

(b)

(d)

(c)

(e)

Figure 1. Molecules used in the present work: (a) SiH4; (b) Si2H6; (c) Si5H12; (d) Si7; (e) Si26. Color codes: white – H, yellow – Si. All electronic structure and MD calculations are performed using the Quantum Espresso (QE) code.95 The core electrons are described by norm-conserving (NC) atomic pseudopotentials. The present QE implementation requires the pseudopotentials of the NC type for computations involving hybrid functionals. Pure functionals can be utilized with the pseudopotentials of different kinds, but we choose NC pseudopotentials to avoid possible differences which arise due to the choice of pseudopotential. Here, we focus only on the differences due to the choice of the functional itself. The present calculations utilize the conditions typically used when performing CPA NA-MD with the neglect of electron-nuclear back-reaction. Under such conditions, the dynamics of electrons is modulated only by the ground-state nuclear dynamics, whereas electronic excitations do not perturb nuclear evolution. Pure density functionals are known to yield reasonably good ground-state properties. Therefore, one can expect the nuclear dynamics computed with pure density functionals to be reasonably close to the one computed with hybrid 8 ACS Paragon Plus Environment

Page 9 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


functionals, especially for condensed-matter system. The closeness further justifies our focus on only the electronic effects. If, however, the electron-nuclear back reaction is not neglected, nuclear dynamics that involves state-specific potential energy surfaces would be different if either pure or hybrid functional were utilized. The valence electrons are described by the Kohn-Sham (KS) orbitals expanded in a planewave basis. The kinetic energy cutoff of 100 Ry is used to achieve the energy convergence of ca. 0.01 Ry for all systems, except for Si26, where cutoff value is reduced to 40 Ry to accelerate the computations. This leads to only a relatively small increase of the energy error, up to 0.015 Ry. The Brillouin zone sampling is performed using only the gamma-point. This is an adequate approach for isolated systems surrounded by a sufficiently large vacuum region. Specifically, a cubic box is used as a unit cell. The size is set to 13 Å for the SiH4 and Si2H6 molecules, 17 Å for the Si5H12 and similarly-sized Si7 molecules, and 20 Å for the Si26 cluster. DFT-MD trajectories for all systems are obtained using the Perdew-Burke-Ernzerhoff (PBE) functional.96,97 The same nuclear trajectory is used to compute NACs with different functionals. This approach helps us single out the effects which arise due to the differences in the electronic wavefunctions and not those associated with differences in nuclear dynamics. This procedure is also consistent with the classical path approximation (CPA)15,54,63 for NA-MD, according to which the nuclear dynamics is predetermined and the electronic degrees of freedom evolve along it. One can envision using a nuclear trajectory predefined by the ground state MD with the PBE functional, but using hybrid functionals to compute the NACs needed to evolve the amplitudes of electronic states. In our calculations, five thousand steps are computed with the 0.5 fs integration time step, resulting in 2.5 ps trajectories. Due to computational limitations, only a part of such



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 36

trajectories can be used in NAC calculations with hybrid functionals. A constant temperature of 300 K is maintained using the Nose-Hoover thermostat. 2.2. Nonadiabatic couplings. The generated MD trajectories are utilized to compute NACs for all pairs of orbitals spanning a several eV energy window around the frontier orbitals. The NACs are computed numerically, using the famous Hammes-Schiffer-Tully formula:  dt  ψ i (t ) ψ j (t + dt ) − ψ i (t + dt ) ψ j (t ) . d ij  t +  = 2 2dt 

(1)

Here, ψ i (t ) is the KS orbital with index i at time t. The dependence of the orbitals on time is implicit, via their implicit dependence on the time-dependent nuclear coordinates, ψ i (R(t )) . The coordinates are dependent on time as computed from the MD trajectories. Within the CPA, the trajectories are obtained using the ground electronic state potential energy surface. The approximation was extensively applied by a number of authors for studying electron and excitation energy transfer in various systems.13,30,73,98–102 We refer readers to the extensive reviews discussing this approach and its applications.18,63,103,104 Due to computational limitations, NACs are always computed at the cost-effective PBE level. The CPA is valid when no significant nuclear reorganization occurs upon electronic excitation. In small silicon QDs, the reorganization energy can vary from 0.1 to almost 1.0 eV, as estimated by various approaches.94,105–107 However, these values correspond to hydrogenated clusters and may be smaller for non-polar bare silicon nanoclusters. In addition, the reorganization energy decreases in large systems, leading to the region where CPA is valid. Under such conditions, the question of how the computed NACs might be affected by Hartree-Fock exchange persists. At the same time, these effects cannot be examined directly in such large systems due to the extremely high computational costs (vide infra). Therefore, the objective of the present work is to evaluate 10 ACS Paragon Plus Environment

Page 11 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


possible differences in NACs computed using pure vs. hybrid functionals for small systems in which benchmarking is practical. Specifically, we consider several functionals available in the QE package. The LDA level is represented by the Perdew-Zunger (PZ) functional,108 whereas the GGA level is represented by the PBE functional,96,97 especially popular for modeling periodic solids. The hybrid HartreeFock/DFT functionals tested comprise the PBE0 (25% of Hartree-Fock exchange)109 and the Heyd-Scuseria-Ernzerhof (HSE06),84,110–113 which is essentially a range-separated version of the PBE0 hybrid. To obtain the statistically-meaningful descriptors, we average the magnitudes of the computed NACs over the trajectory snapshots: d i, j =

1 N

N −1

∑

d i , j (t k ) ,

(2)

k =0

where N is the number of the MD configurations included in the averaging. This number varies, because of the extremely demanding computations involving hybrid functionals in large systems, like Si26. Specifically, all 5000 configurations are utilized in the calculations involving PZ and PBE functionals, whereas we use 1000 – 1500 configurations if the calculations are performed with hybrid functionals.

3. Results and discussion The NACs computed for several considered systems are shown in Figure 2. We illustrate the results for the smallest (SiH4 molecule), the intermediate (Si7 cluster), and the largest (Si26 quantum dot) systems studied in the present work. The results for the pure (PBE) and two hybrid (PBE0 and HSE06) functionals are presented. Note also that the HSE06 functional incorporates the range-separation feature, which is especially important for minimizing the delocalization 11 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 36

error in larger systems.88 One can observe two main trends. First, NACs computed using the PBE functional are notably different from those obtained from the calculations employing the hybrid schemes. On average, the inclusion of the Hartree-Fock exchange leads to decreased couplings. Second, the observed difference is diminished in small molecular systems and increases in larger clusters. Both PBE0 and HSE06 functionals lead to comparable NAC magnitudes. Small difference develops as the system’s size increases. Yet, this difference is negligible in comparison to the difference due to the presence of the exact exchange. Thus, the range separation of the added Hartree-Fock exchange does not appear to be as critical as its nonnegligible presence. The difference between NACs computed using PBE0 and HSE06 functionals may develop as the system’s size increases.

(a)

(b)

(d)

(e)

(c)

(f)


Page 13 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(g)

(h)

(i)

Figure 2. Non-adiabatic couplings computed with three density functionals: PBE (a, d, g), PBE0 (b, e, h), HSE06 (c, f, i) for three test systems: SiH4 (a, b, c), Si7 (d, e, f), and QD-Si26 (g, h, i). The visualization in Figure 2 is a convenient way to understand main qualitative differences. The comprehensive quantitative characterization of the computed couplings is presented in Table 1. In addition, the couplings are sorted into three groups, depending on the type of orbitals in the pair: a) occupied-occupied (OO), these couplings are relevant to hole relaxation dynamics; b) occupied-unoccupied (OU), these couplings are related to nonradiative electron-hole recombination rates; c) unoccupied-unoccupied (UU), these couplings are responsible for photoexcited electron relaxation dynamics. For small silicon hydride molecules, the average couplings in the OO or UU blocks are notably larger than in the OU block. In the silicon-only systems, the couplings in all blocks are comparable to each other. Table

1.

Average

NAC

(meV)

for

the

occupied-occupied/occupied-

unoccupied/unoccupied-unoccupied orbital pairs. SiH4

Si2H6

Si5H12

Si7

Si26

PZ

40.8/1.0/22.7

21.5/1.0/22.3

43.2/0.5/38.8

4.5/3.8/5.0 9.1/8.3/7.8

PBE

40.5/0.9/25.3

22.2/0.8/21.7

45.0/0.5/40.2

5.6/6.3/6.1 13.0/19.1/12.1



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

PBE0

50.0/0.6/32.0

21.9/0.4/23.7

23.6/0.2/17.2

0.4/0.3/0.7 1.5/2.6/1.5

HSE06

46.8/0.7/30.6

23.6/0.5/22.7

18.6/0.2/15.1

0.7/1.0/0.6 1.3/1.9/1.0

Page 14 of 36

Two trends inferred from Figure 2 can be further analyzed on the basis of data summarized in Table 1. First, one can clearly see that couplings in the OO and UU blocks are notably larger in silicon hydrides than they are in bare silicon clusters. This effect is not necessarily related to the mere size of the molecular system, as can be seen from the comparison of the couplings in Si2H6 to the NACs in the other two (smaller and larger) hydrides. More likely, the presence of polar SiH bonds with high vibrational frequencies induces larger couplings in the hydrides, as can be explained by the bond polarity and vibrational frequencies. Large polarity of the bond implies large dissimilarity of the involved orbitals, leading to large response of their overlap to external perturbations (e.g. nuclear displacements). In contrast, bare Si clusters are composed of alike atoms, leading to a reduced bond polarization. The MOs in these systems are well delocalized over the entire cluster and are less sensitive to nuclear motion. High Si-H vibrational frequencies increase the rate with which MO overlaps can change in time. On the contrary, Si-Si vibrations in pure Si clusters are slower, leading to smaller NAC magnitudes. Both factors can be understood on the basis of Eq. 1. The second trend observed is the reduced magnitude of the NACs computed with the hybrid functionals. As can be seen from Table 1, the NACs computed for SiH4 and Si2H6 with PBE or PZ functionals are on the same order of magnitude as those computed with PBE0 or HSE06. As one moves toward larger systems, the difference (or, rather, the ratio of) couplings increases: in Si5H12 the values vary by a factor of 2, whereas in Si7 or Si26, they differ by up to an order of magnitude. The qualitative explanation of this effect can be obtained by understanding the


Page 15 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


significance of the exact exchange in the hybrid DFT schemes. By the construction, hybrid functionals substitute a fraction of the pure (based only on the electron density) exchange energy with the exact exchange energy from the Hartree-Fock theory (hence requiring explicit orbitals). In the Hartree-Fock theory, the exchange energy appears from the Slater determinant construction that reinforces the Pauli repulsion of electrons. Thus, the exchange is responsible for minimizing the clashes of electrons existing in different states (orbitals). In other words, the exchange contribution precludes an excessive delocalization of orbitals (and hence affecting the shapes of MOs). Localized orbitals have smaller chances to develop large NACs, as also suggested by Eq. 1. Consequently, the exact exchange contribution minimizes NACs. Having said that, we can recall that the native DFT exchange functionals are short-ranged. As the system’s size increases, this shortcoming becomes more and more critical, allowing distant electrons to clash more than they are supposed to. As a consequence, the overlaps of the orbitals at adjacent times, ψ i (t ) ψ j (t + dt ) , and, therefore, NACs between the involved orbitals increase on average. On the contrary, the Hartree-Fock exchange is a long-ranged term. Thus, it would enforce even the electrons separated by a relatively large distance to obey the Pauli exclusion principle. As such, the coupling of the orbitals decreases. In view of the above discussion, it is now clear why the differences between pure and hybrid functionals start manifesting themselves only in sufficiently large systems. Both SiH4 and Si2H6 are likely be of the size suitable for proper exchange description by the pure PBE functional. Not so for Si5H12, Si7, and especially not for Si26 molecules. The size dependence of NACs found here is also consistent with the trends reported recently by the Isborn group.88 It was shown how the delocalization error present in pure DFT functionals increases as the system size grows, affecting the computed electronic properties. Their results



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 36

also suggest that global hybrids may still show a notable delocalization error, whereas rangeseparated hybrids may be more favorable. Our calculations show that the difference in NACs computed using the global hybrid, PBE0, and the range-separated one, HSE06, are not that drastic. The distinction may increase as one goes to even larger systems, especially if the electronic states are well-delocalized or extended. In that region, the full exact exchange added for large inter-electron separations in long-range corrected hybrids may become important for ensuring the Pauli exclusion principle. Based on the present results, we expect that such functionals would lead to a further decrease of the NAC magnitudes, just as hybrids of both types decrease the NACs in comparison to pure functionals, when applied to the small systems considered here. The differences in computed NACs can also be understood from the simple mathematical point of view. On the one hand, NACs are proportional to the derivative coupling vectors,

d ij ~ ψ i

d ψ j , ∀α . On the other hand, using Hellman-Feynman theorem, one can show that dRα

the derivative couplings are inversely proportional to the energy gap between the adiabatic states:

d ij ~ ψ i

d ψj ~ dRα

ψi

dH ψj dRα , ∀α . Ei − E j

(3)

Pure density functionals are notorious for their underestimation of the electronic gaps, Ei − E j , for the same physical reason as discussed above. As a result, the underestimated gaps can lead to the overestimated NACs. The situation is likely to be more dramatic, the closer the PBE (or PZ) gap is to zero. To verify this line of thought, we compute the projected density of states (pDOS) for the considered systems (Fig. 3).


Page 17 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 3. Projected densities of states in SiH4 (a, b, c), Si7 (d, e, f), and Si26 (g, h, i) computed using PBE (a, d, g), PBE0 (b, e, h), and HSE06 (c, f, i) functionals. Peaks are broadened by



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 36

convolving the computed energy levels with the Gaussian distribution. Arbitrary units for the y axis are chosen. As it is expected, the HOMO-LUMO gap decreases as the size of the system increases. For SiH4 (Fig. 3 a-c), the computed gaps vary in the 7.7-9.4 eV range (7.7 with PBE, 9.4 - PBE0, and 8.8 eV – HSE06). The PBE value is relatively close to the PBE0 and HSE06 gaps. Consequently, there are no notable differences in NAC values. As we move toward the larger (and nonhydrogenated) Si7 cluster, the differences in the HOMO-LUMO gap become more notable. The PBE calculations predict only a 0.3 eV gap, whereas PBE0 and HSE06 predict significantly larger gaps of 1.4 eV and 0.7 eV, respectively. Finally, the Si26 cluster has a vanishing gap of 0.2 eV if computed with the PBE functional, whereas the PBE0 and HSE06 functionals open it up to 1.1 eV and 0.5 eV, respectively. The HOMO-LUMO gaps computed for all systems are summarized in Table 2. The existing theoretical and experimental reference values are also given. We should note, that the gap values for small Si clusters may range significantly, depending on the methodology used and the details of the molecular structure. Although, various calculations show gaps in Si clusters up to 3.5 eV,

114

the existing scanning tunneling spectroscopy

experiments115,116 suggest that only silicon clusters smaller than 15 Å in diameter show nonvanishing gaps up to 0.45 eV. Such small gaps originate from the presence of surface defects and dangling bonds. The gaps in our PBE/PZ calculations vary around 0.2 eV. They are increased up to 1.4 eV and only up to 0.5-0.7 eV when the PBE0 and HSE06 functionals are employed, respectively. The computed gaps fall reasonably within the ranges reported by previous studies, which demonstrates the applicability of the utilized methods for describing the electronic structure of Si clusters and Si hydrides. Bare silicon clusters represent an interesting type of system – they represent clusters with very small energy gaps, yet not metallic in nature. For such 18 ACS Paragon Plus Environment

Page 19 of 36

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


systems, the differences in electronic structure properties (gaps, NACs, etc.) computed with pure and hybrid functionals are expected to be more notable than in the systems with large HOMOLUMO gaps. This expectation is indeed observed in our computations (Fig. 2, Table 1).

Table 2. HOMO-LUMO gaps (eV) for the PBE-optimized structures of all considered systems. SiH4

Si2H6

Si5H12

Si7

Si26

PZ

7.8

6.5

5.6

0.2

0.2

PBE

7.4

6.4

5.7

0.3

0.2

PBE0

9.4

7.9

7.5

1.4

1.1

HSE06

8.8

7.3

6.8

0.7

0.5

Reference

7.98.3117,118

6.767.17120

ca. 5.7119

Expt.

Dependence of Nonadiabatic Couplings with KohnâSham Orbitals on

Recommend Documents

Dependence of Nonadiabatic Couplings with KohnâSham Orbitals on