Scaling of the self-energy correction to the HOMO-LUMO gap with magnesium cluster size and its potential for extrapolating to larger magnesium clusters Free

Accurate computations of the electronic structures are critical for screening and design of materials for various applications. Density function theory (DFT) with the local density approximation (LDA) or the generalized gradient approximation has been successful in predicting the ground-state properties of materials.¹ DFT, however, is inapplicable for calculating accurate electronic energy levels, as the simplification boosting the computational efficiency neglects some important features of the self-energy term ∑ of electrons. Consequently, DFT always underestimates the bandgaps for systems with energy gaps across the Fermi levels, such as metallic clusters. As an accurate approach, the GW quasiparticle (QP) method has been derived based on the many-body Green's function theory (MBGFT).² The method preserves the non-local and dynamic features of ∑ of electrons and has been proven to be applicable for calculations of electronic structures of different physical systems, including bulk materials, molecules, surfaces, and clusters.^3–5 Numerous predictions of electronic structures of solids and organic molecules by the GW method agree well with experiments.^6–11 More recently, Marom and coworkers successfully computed ionization potentials and electron affinities of organic molecules.^12,13

In spite of the tremendous success of the GW method in computational materials, it has two drawbacks: the discrepancies of results for the same material by different researchers and the enormous amount of computational cost.^14–16 The large deviations of published results were attributed to the basis set convergence of GW calculations, the accuracy of the pseudopotential approximation, and the various choices for the starting mean-field calculations.¹⁷ In a recent benchmark investigation, researchers found that partially or non-self-consistent GW methods are often more accurate than fully self-consistent GW, which was attributed to error cancellation, and a general solution for complex systems such as nanostructures had yet to be found.¹⁸ 

The computational cost of a typical GW algorithm is generally O(N⁴) with respect to system size N, with a prefactor.¹⁹ To resolve the vast computational cost, researchers have proposed to reduce the prefactor by avoiding summation over empty states in the polarizability and a low-rank approximation of the dielectric matrix.²⁰ Foerster and coworkers reduced the computational complex to O(N³) by using a cubic-scaling GW algorithm in a Gaussian basis.²¹ Another study devised a stochastic GW algorithm, which could scale linearly with the system size. The method has been applied to a silicon cluster with a thousand atoms, and application to more complex systems is underway.^22,23

Observing that the density of states (DOS) of materials resembles that of the free electron gas at high energies, Gao and coworkers simplified the summation over high-energy conduction in GW calculations to be a numerical integration on a sparse energy grid.²⁴ The method enables fully converged GW calculations for large systems feasible, eliminating the need for performing expensive and tedious convergence tests, where all conduction bands in GW calculations can be included effectively at greatly reduced computational costs. They successfully demonstrated their method for oxides of ZnO and MgO supercells. In another recent investigation, Wilhelm and coworkers presented a full-frequency GW algorithm in a Gaussian-type basis, with a scaling of O(N²) to O(N³).²⁵ The method adopts the GW space-time method with sparse linear algebra and minimax grids for imaginary time and frequency.^26–28 

Yet almost all the above methods and computations require massively parallel execution on state-of-the-art supercomputers.^{17,18,23–26} While the development of new algorithms within the GW method to reduce the computational cost of the GW method is always warranted, new methods of computation with high accuracy of GW and low cost are also desired. In this paper, we propose a methodology for calculating energies of HOMOs and LUMOs for clusters with the accuracy of the GW method, while with the computational feasibility of DFT. In the computation, the energy of a single electronic state is calculated as $EDFT+Δ$⁠, where $EDFT$ is the energy of this state obtained by DFT and $Δ=EGW−EDFT$ is usually referred to as “scissors correction.” The algebra does not make any substantial progress, if one still needs the GW implementation to determine $EGW$ and $Δ$⁠. However, if there is a solid approximation $Δ$′ for $Δ, EDFT$ alone can be used to determine the electronic energies accurately. The energy differences between $EGW$ and $EDFT$ in some bulk materials have been investigated,²⁹ and it is found that $Δ$ is energy dependent. It could be much more complicated for clusters with different sizes, since $Δ$ could also be size dependent. The investigation finds that the energy difference of the scissors correction is independent of energy, and thus, the HOMOs and LUMOs can be efficiently calculated using DFT only.

In the GW implementation, the energy of an electronic state is obtained by solving the quasiparticle (QP) equation

where $T$ is the kinetic energy operator, $Vext$ the external potential, $VH$ the Hartree potential, and $Ei$ and $ϕi$ the energy and wavefunction of the ith QP. $Σxc(r,r′;Ei)$ is the exchange-correlation self-energy operator, which includes the non-local and dynamic effects of the electronic interaction. The QP equation is solved based on the results of the density functional theory (DFT)

where $εi$ and $φi$ are the eigenvalue and eigenfunction of the ith Kohn-Sham (KS) particle, respectively, and $Vxc(r)$ the exchange-correlation potential. With the assumption that the KS eigenfunctions $φi$ agree well with the QP wavefunctions $ϕi$ in most cases,^30,31 QP energies are usually solved with perturbative method to the first order

According to Hedin's equations,³² $Σxc=iGWΓ$ (with $ℏ=1$⁠), where Г is the vertex function and $G$ is the one-particle Green's function

The coefficient $ηn$ is $+1$ for unoccupied states and $−1$ for occupied states. W is the screened Coulomb interaction which can be written as $W=V+VΠV$⁠, where $V(r,r′)$ is the Coulomb interaction, and $Π(r,r′;E)$ is the reducible polarizability which represents the linear response of the electronic density to the external perturbative potential and can be obtained by solving the equation $Π=P+PVΠ$⁠, where $P=−iGGΓ$ is the irreducible polarizability. Both $Σxc$ and $Π$ (or W) include the vertex function Г. It has been shown that a consistent choice of Г is necessary for the QP calculation.³³ The present calculation goes beyond the standard RPA based GW approximation by including vertex corrections to both the self-energy $Σxc(r,r′;Ei)$ and the reduced polarizability $Π$ in terms of the time-dependent local density approximation (TDLDA) exchange correlation kernel. The relevant matrix elements of the vertex part of the self-energy are essentially obtained from the solution of the Cassida equations.³⁴ More details can be found in Ref. 35.

DFT calculations for the ground state of all magnesium clusters Mg_n (n = 3–22, 25, 30, 35, and 40) are performed using the SIESTA¹ method and code with the local density approximation (LDA) as the exchange and correlation functional. The core electrons [1s²2s²2p⁶] of Mg are replaced by the nonlocal norm-conserving pseudopotential based on the Troullier-Martins scheme,³⁶ and the Ceperley and Adler³⁷ correlation as parametrized by Perdew and Zunger² is used. A double-ζ polarization (DZP) basis set of numerical atomic orbitals is used for the valence electrons of Mg, with cutoff radii of 10 a.u. for both s and p orbitals.

Detailed calculations for the matrix elements of Eqs. (2) and (3) involving potentials can be found in Ref. 1, which are performed on a real-space grid. Briefly, the eigenfunctions are expressed in the atomic basis set, and the electron density is computed accordingly. When the electron density is found, one solves Poisson's equation to obtain the three terms consisting of total grid potential V(r). The exchange integrals $∫dr∫dr′φi(r)φj(r)V(r,r′)φk(r′)φl(r′)$ are then evaluated with the multigrid method.⁴ 

The GW$Γ$ scheme is implemented in this paper due to its better performance over the GWA scheme (⁠$Γ=1$⁠) for simulations of Mg clusters.³⁸ The structures of Mg clusters are optimized by simulated annealing using molecular dynamics with an exponential cooling schedule, followed by the conjugated gradient algorithm with the maximum force tolerance of 0.01 eV/Å. The optimized structures of Mg clusters are presented in Fig. 1. All integrals are evaluated on a uniform grid in real space with a grid spacing of 0.6 a.u., which has been tested to give QP energies with an accuracy of 0.1 eV. The convergence of the QP calculation usually requires a large number of unoccupied states for the evaluation of the polarizability. The convergence of the correlation part $⟨φi|Σc|φi⟩$ can be accelerated by truncating the summation at some point and evaluating the reminder within the COHSEX (Coulomb hole + screened exchange) approximation.³⁵ The COHSEX approximation is essentially a static limit of $⟨φi|Σc|φi⟩$⁠, which assumes that screening is instantaneous, and the polarizability and the potential W₀ become proportional to a δ-function in time and constant in energy. The integral over energy can then be replaced by a sum over poles of the single-particle Green's function. In the implementation, one can recover the COHSEX approximation by imposing $ωs≫|E−εi|$ and calculation of the reminder.³⁵ 

Note that wavefunctions $ϕi$ of a finite system are real functions, while the corresponding energies $Ei$ are complex. The imaginary parts of $Ei$ arising from the dynamic character of the self-energy operator are related to the scattering rates and finite lifetimes of electronic states.³⁴ In this sense, the QP equation [Eq. (1)] is essentially a complex equation. With the complex QP energies written explicitly as $Ei−iηiγi$⁠, Eq. (3) is now expressed as

which is equivalent to an equation set composed of two real equations

It should be emphasized that each electronic state corresponds to a single equation set. However, all equation sets of a system with multiple electronic states are correlated, since the term $Σxc(Ei−iηiγi)$ depends explicitly on all energies $En$ of electronic states. Therefore, all energies $En$ have to be determined simultaneously with numerical iteration. This approach has been used to calculate energies and lifetimes of electronic states in Mg and Si clusters.^38,39

In the QP calculations, a ready starting point is the approximation for G

From $G0$⁠, one can evaluate $Σxc$ as $iG0WΓ$ and thus solve Eq. (1) for QP energies $En$⁠. We have demonstrated that it is crucial to calculate iteratively until the convergence of $Ei$ to make the GW implementations justified physically and numerically.^38,39 Blase and co-workers' recent investigation also demonstrated that self-consistent schemes provide a good agreement with reference calculations.⁴⁰ All GW calculations in this paper are performed self-consistently.

Since the convergence could be an important issue for the GW implementation^17,41, both the double-zeta polarization (DZP) basis set and the quadruple/quintuple-zeta polarization (4Z4P) basis set are tested to evaluate the influence of the completeness of the basis set for the GW simulations of Mg clusters in the investigation. The calculated electronic affinities obtained with DZP and 4Z4P basis sets for Mg10 are 1.49 eV and 1.45 eV, respectively, or the difference is within 0.05 eV. The calculated ionization potential with DZP for Mg10 is 5.47 eV. Those results agree well with experiments: The experimental value of the electron affinity Mg10 was found to be 1.6 eV, and the ionization potential was 5.43 eV.⁴² This means that both DZP and 4Z4P basis sets are the complete basis set for the magnesium clusters. As presented in Sec. III B, the computed lowest unoccupied molecular orbitals (LUMOs) of $Mgn−$ clusters are in agreement with the first photoelectron spectra (PES) peaks of $Mgn−$ clusters, which shall further validate the applicability of the double-ζ polarization (DZP) basis set.

The highest occupied molecular orbitals (HOMOs) and the lowest unoccupied molecular orbitals (LUMOs) of Mg_n (n = 3…22, 25, 30, 35, and 40) obtained by the LDA and GW methods are given in Table I, together with the scissors operators $Δ$ calculated as the energy differences. Those results in Table I can be comprehensively compared to the measured the photoelectron spectra (PES) of $Mgn−$ clusters.⁴² 

For example, for Mg₅, Mg₁₀, Mg₁₅, Mg₂₀, Mg₂₅, and Mg₃₀, the GW computed LUMOs/experimental electron affinities are 0.92/1.15 eV, 1.49/1.6 eV, 2.11/2.05 eV, 2.06/1.95 eV, 2.53/2.4 eV, and 2.55/2.45 eV, respectively. Thus, the LUMOs calculated by the GW method in this paper are in agreement with the first PES peaks of $Mgn−$ clusters, which can be used as the reference for the LUMOs of corresponding neutral Mg clusters. This further demonstrates the convergence of the computations by using the double-ζ polarization (DZP) basis set for the Mg clusters.

The comparison between HOMOs and LUMOs of Mg_n(n = 3…22) is also plotted in Fig. 2(a). In general, the HOMO increases while the LUMO decreases with the increasing particle size, which is consistent with previous published results obtained by DFT.^43,44 It can also been seen that the GW correction consistently shifts up the LUMOs and shifts down the HOMOs of the simulated Mg clusters and therefore widens the bandgaps obtained by DFT.

According to the electronic shell model, delocalized valence electrons in a metallic cluster will fill the cluster in a similar way as in a single atom.³⁴ Such an atom-like shell filling effect has been demonstrated to be the dominant factor on the size abundance of small magnesium clusters produced in superfluid helium droplets.⁴⁵ The effect also manifests itself in the DFT and GW results shown in Fig. 2(a), where the LUMO of Mg₄ is notably low when compared with Mg₅, indicating the electronic shell closure in Mg₄ (1p shell). HOMOs in Fig. 2(a) undergo a large jump from Mg₉ to Mg₁₀ and from Mg₁₉ to Mg₂₀, which means that Mg₁₀ and Mg₂₀ have stable electronic states also attributed to the electronic shell closure (2s and 2p shells). The differences between the DFT results and GW results are calculated as the scissors operator $Δ$ as plotted in Fig. 2(b). It is notable to see that $Δ$ is rather insensitive to the electronic shell effect, which has been smoothed out when calculating the energy difference. This indicates that the scissors operator for HOMOs and LUMOs can be approximated as a function of the cluster size only, with the complicated electronic shell effect being neglected.

To quantify the dependence of the scissors operator on the cluster size, $|Δ|$ of HOMOs and LUMOs of Mg_n (n = 3–22, 25, 30, and 35, 40) is plotted versus n in a log-log style in Fig. 3. We find that both datasets demonstrate sound linearities, as indicated by the two lines in Fig. 3. This implies that the size-dependence of $|Δ|$ can be described by a simple power function

where the coefficients a and b determined through the linear fitting are 2.79 and −0.34 for $|Δ|$ of HOMOs and 2.38 and −0.24 for $|Δ|$ of LUMOs. With such an empirical correction, we can predict the HOMOs and LUMOs of larger metallic clusters with the GW accuracy using a DFT computation. Note that $Δ$ of HOMOs and LUMOs should vanish for infinite large n, or the bulk magnesium will be predicted to have a non-vanishing bandgap after the scissors operators. This requirement is satisfied implicitly during the linear fitting, as the two power coefficients obtained above are negative.

A more interesting issue about the scissors operators is whether $Δ$ of HOMOs (LUMOs) estimated as above can be applied to all occupied (unoccupied) states, since most scissors operators in literatures are introduced as a rigid energy shift over the DFT energies for all unoccupied states, instead for the HOMO and LUMO only, to reproduce bandgaps measured experimentally. In Fig. 4, $Δ$ of occupied and unoccupied states around the Fermi levels is plotted versus the DFT energies for some Mg clusters (Mg₂₀, Mg₃₀, and Mg₄₀). It is found that for each Mg cluster, the scissors operators vary little for occupied (unoccupied) states around the Fermi level. This justifies the application of a rigid energy shift over the DFT energies for occupied (unoccupied) states near the Fermi level in metallic clusters. It should be pointed out that the approximation is validated only for states near the Fermi level, with which most electronic and optical properties of metallic clusters are associated. For states far below or above the Fermi level, a rigid energy shift may deviate from the accurate GW results significantly and thus shall be inappropriate.

In summary, we have simulated electronic structures of magnesium clusters Mg_n (n = 3–22, 25, 30, 35, and 40) using the density functional theory (DFT) and the quasiparticle GW method. The scissors operators (GW corrections over DFT energies) are found to be size-dependent only with a negligible complicated electronic shell effect and thus are fitted by simple power functions of the cluster size n. The highest occupied molecular orbitals (HOMOs) and the lowest unoccupied molecular orbitals (LUMOs) of larger clusters can be accurately calculated using DFT computation and the scissors correction. It is also demonstrated that the scissors operators are appropriate for occupied and unoccupied states around the Fermi level.

This project was sponsored by the National Science Foundation (Grant No. CBET-0830098) and Hunan provincial grant.