Linear response approach to collective electronic excitations of solids and surfaces

https://doi.org/10.1016/j.cpc.2008.10.019Get rights and content

Abstract

We have developed a parallel computer program for the study of dynamic response of periodic systems. It computes the linear response of an interacting many-electron system from its ground-state electronic structures, which are obtained from ab initio band structure calculations in the plane-wave and pseudopotential scheme. As test examples, we applied this program to calculate the linear response of bulk aluminum and a beryllium monolayer. The excitation spectra show prominent plasmon resonances, which compare well with the available data and previous calculations. For surfaces or thin films, we found that removing periodicity perpendicular to the surface gives a more reliable description of the low-energy excitation spectra, especially in the long-wavelength limit.

Program summary

Program title: Dresponse

Catalogue identifier: AECK_v1_0

Program summary URL: http://cpc.cs.qub.ac.uk/summaries/AECK_v1_0.html

Program obtainable from: CPC Program Library, Queen's University, Belfast, N. Ireland

Licensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.html

No. of lines in distributed program, including test data, etc.: 49 098

No. of bytes in distributed program, including test data, etc.: 11 836 088

Distribution format: tar.gz

Programming language: Fortran 90/MPI

Computer: Any architecture with a Fortran 90 compiler

Operating system: Any

Has the code been vectorized or parallelized?: Yes

RAM: 50 MB–2 GB per processor depending on system size

Classification: 7.3

External routines: BLAS (http://www.netlib.org/blas/), Lapack (http://www.netlib.org/lapack/), MPI (http://www-unix.mcs.anl.gov/mpi/), abinit (for ground-state calculations, http://www.abinit.org/)

Nature of problem: The dynamic response of bulk and surface systems. It is usually dominated by collective electronic excitations (plasmons) at low-energy range.

Solution method: The ground-state wavefunctions are obtained from ab initio density-functional calculation in the planewave and pseudopotential scheme [1]. The linear response theory combined with the time-dependent density functional theory is implemented for the investigation of excitation properties.

Restrictions: The present version only handles 3D and 2D periodic systems.

Unusual features: A mixing reciprocal/real-space basis is implemented in surface calculations in order to remove the intercell coupling in solving the Dyson equation. This treatment provides more reliable results, especially in the long-wavelength limit.

Running time: The example included in the distribution takes a few minutes to complete.

References:

[1] G. Onida, L. Reining, A. Rubio, Rev. Mod. Phys. 74 (2002) 601.

Introduction

Linear response theory combined with time-dependent local density approximation (LR-TDLDA) [1], [2], [3], [4] provides one of the most accurate approaches to the study of dynamic response of interacting many-electron systems. It has widely been applied to optical properties of finite systems like atoms [1], [2], molecules [5] and small metal clusters [6], [7], optical and electron energy-loss spectroscopy (EELS) of solids and their surfaces [8], [9]. For bulk metals [10] and their surfaces [11], this approach has been very successful in the description of collective electronic excitations, the plasmon, which usually dominates the dynamic response of a many-electron system.

The electron-gas or the jellium model offers a simple approach to excitations of valence electrons in most free-electron-like metals, where TDLDA has been able to qualitatively reproduce the frequencies and dispersions of bulk and surface plasmons [12], [13], [14]. For noble and transition metals, the jellium approximation has been problematic due to the tight-binding character of the d electrons. Even for simple metals, the jellium model results have been found to deviate quantitatively from experiment. For example, it cannot describe the positive linewidth dispersion of potassium plasmon [15], and the finite lifetime of magnesium surface plasmon in the long-wavelength limit [16], [17]. Such differences can be attributed to the band structure effects, which were neglected in the jellium approximation.

The discrepancies between the jellium model and experiments could be dramatically reduced in atomistic modeling. Quong and Eguiluz first took a band-structure approach to LR-TDLDA using the planewave-pseudopotential (PWPP) method and found that the plasmon dispersions of aluminum and sodium were in better agreement with experiment than the jellium model [18]. Later, Aryasetiawan et al. implemented a linear response scheme using the linear muffin-tin orbital (LMTO) basis and affirmed the measured negative energy dispersion of Cs plasmon [19], [20]. More recently, Ku and Eguiluz developed a TDLDA scheme based on the accurate full-potential-linearized-augmented-planewave (LAPW) method [15]. The positive linewidth dispersion of potassium plasmon was found to result from the decay channel provided by the unoccupied d bands, which are more than 4 eV above the Fermi level. The first study of surface excitations with ab initio band structure were carried out by Silkin and collaborators, who calculated the dynamic response of metal monolayers [21] and lifetimes of surface plasmons [17]. So far, dynamic responses of a number of solids and surfaces have been reported extensively [22], [23], [24], [25], [26], [27], [28], [29], [30], [31].

TDLDA has been implemented into many softwares, such as abinit [32], cpmd [33], octopus [34], etc. Most of them adopt the Casida's approach [35], [36], which is only valid for finite systems. To the best of our knowledge, the only open-source TDLDA software applicable for extended systems is dp [37], where the computation is carried out in reciprocal space. For surfaces and thin films, the present implementation in reciprocal space [21], [37] overestimates the long-range Coulomb interaction between neighboring unit cells normal to the surface. This artificial effect becomes prominent in the long-wavelength limit [21].

Here, we present a LR-TDLDA program for the calculation of dynamic response of periodic systems. This program was mainly designed to provide efficient numerical approach to surface plasmons in low-dimensional structures. For this purpose, we adopt a single-slab geometry in real space. We found that our real-space scheme, which removes the long-range Coulomb interaction between neighboring unit cells, gives more reliable results and is computationally more efficient, especially in the long-wavelength limit. The program is fully parallelized using message passing interface (MPI). Currently, the input band structure and wavefunctions are taken from abinit [38], [39], a density functional theory (DFT) [40], [41] package.

The rest of this paper is organized as follows. In Section 2, the theoretical method of LR-TDLDA is briefly formulated. The details of the code implementation and parallelization are given in Section 3. The structure of the program is outlined in Section 4. In Section 5, we present two examples of applications, bulk aluminum and a beryllium monolayer, where comparison with previous calculations and experiments are given. The memory requirement and parallel scaling in computations are discussed. It is followed by the discussion of the validity of TDLDA description of plasmons. Finally, a summary is given in Section 6. Atomic units e2=m==1 are used in this paper unless stated otherwise.

Section snippets

Linear response theory and TDLDA

The central quantity in LR-TDLDA is the density response function, χ(r,r,ω), which is defined by [42]δn(r,ω)=drχ(r,r,ω)Vext(r,ω), where δn is the electron density induced by external potential Vext. This response function embodies many excitation properties of an interacting electron system in the linear response approximation.

Within either time-dependent density functional theory (TDDFT) or the many-body perturbation theory, χ of an interacting system can be obtained from its

Symmetry operations on Bloch states

Ground state DFT calculations are usually performed by a small set of k-points in the irreducible Brillouin zone (IBZ). However evaluation of the response function χ0 in Eq. (7) requires summation of all k-points in the whole BZ. To minimize the computation and data storage, the Bloch states outside the IBZ can be generated from those inside the IBZ through symmetry operations [44]Ψn,R˜1k(r)=eikτRΨn,k(R˜r+τR), where R˜ is a rotation operator and τR is the translation vector associated with

Structure of the program and input files

Fortran source code files of the program are listed in this section.

  • ab.f90 is the main program which controls each step of the calculation.

  • rdabinit.f90 reads the wavefunction file output by abinit.

  • pwbasis.f90 initializes the planewave basis in the linear response calculation. Note that this basis can be different from the one used in the ground state calculation.

  • kpoints.f90 finds all k-points in the whole BZ from those in the IBZ using the symmetry operations.

  • chi0.f90 constructs χG,G0(q,ω).

Test examples

We present two test examples, bulk aluminum and a beryllium monolayer, for which comparisons with early calculations or experiments are possible. For the bulk calculation, our method fully reproduces the earlier theoretical results and measured data. In the surface case, our G-space results are also in agreement with early calculations [21]. In addition, it is shown that our real-space approach gives a more accurate description of the low-energy excitations, especially in the long wavelength

Summary

We have developed a parallel computer program for the calculation of the dynamic responses of solids and surfaces, based on the ab initio band structure and wavefunctions. Test calculations of excitation spectra for bulk aluminum and a beryllium monolayer are in close agreement with previous calculations. In particular, for the beryllium monolayer, we have demonstrated that the real-space implementation is able to remove the intercell coupling and provides more reliable description of

Acknowledgements

We thank V.M. Silkin for helpful discussions. This work has been supported by the PhotoNano frame project, the material consortium-ATOMICS, funded by the Foundation for the Strategic Research (SSF), the Institutional Grant Program of STINT, Sweden, and the 100-Talent Program of Chinese Academy of Sciences.

References (66)

  • A.G. Eguiluz et al.

    Nucl. Instrum. Methods B

    (1995)
  • V.M. Silkin et al.

    Vaccum

    (2006)
  • X. Gonze

    Comput. Mater. Sci.

    (2002)
  • S. Goedecker

    Comput. Phys. Comm.

    (1993)
  • Z. Yuan et al.

    Surf. Sci.

    (2008)
  • A. Zangwill et al.

    Phys. Rev. A

    (1980)
  • M.J. Stott et al.

    Phys. Rev. A

    (1980)
  • M.A.L. Marques

    Time-Dependent Density Functional Theory

    (2006)
  • G. Onida et al.

    Rev. Mod. Phys.

    (2002)
  • Z.H. Levine et al.

    Phys. Rev. A

    (1984)
  • W. Ekardt

    Phys. Rev. B

    (1985)
  • D.E. Beck

    Phys. Rev. B

    (1984)
  • S. Botti et al.

    Rep. Prog. Phys.

    (2007)
  • J.M. Pitarke et al.

    Rep. Prog. Phys.

    (2007)
  • A. Liebsch

    Electronic Excitations at Metal Surfaces

    (1997)
  • K. Sturm

    Adv. Phys.

    (1982)
  • K.D. Tsuei et al.

    Phys. Rev. Lett.

    (1989)
  • K. Tsuei et al.

    Phys. Rev. Lett.

    (1991)
  • W. Ku et al.

    Phys. Rev. Lett.

    (1999)
  • H. Ishida et al.

    Phys. Rev. B

    (1996)
  • V.M. Silkin et al.

    Phys. Rev. Lett.

    (2004)
  • A.A. Quong et al.

    Phys. Rev. Lett.

    (1993)
  • F. Aryasetiawan et al.

    Phys. Rev. Lett.

    (1994)
  • F. Aryasetiawan et al.

    Phys. Rev. B

    (1994)
  • A. Bergara et al.

    Phys. Rev. B

    (2003)
  • N.E. Maddocks et al.

    Phys. Rev. B

    (1994)
  • I. Campillo et al.

    Phys. Rev. B

    (1999)
  • M.A. Cazalilla et al.

    Phys. Rev. B

    (2000)
  • A.G. Marinopoulos et al.

    Phys. Rev. Lett.

    (2002)
  • N. Vast et al.

    Phys. Rev. Lett.

    (2002)
  • W. Ku et al.

    Phys. Rev. Lett.

    (2002)
  • W.D. Schone et al.

    Phys. Rev. B

    (2003)
  • I.G. Gurtubay et al.

    Phys. Rev. B

    (2005)
  • Cited by (27)

    View all citing articles on Scopus

    This paper and its associated computer program are available via the Computer Physics Communications homepage on ScienceDirect (http://www.sciencedirect.com/science/journal/00104655).

    View full text