Lausanne 2017Version 4 out and rockingYambo-pyFLASH-IT

Main menu

HomeNewsPeopleDownloadRun the codeInput fileTutorials
Overview GW Lifetimes SiH4 Fantastic dimensions Hydrogen chain Electron Phonon KERR effect Surface spectroscopy GaSb Parallel Developing Yambo
DocumentationPublicationsEventsContactsRobots Forum

The Wikipedia Page of Yambo Yambo@Wiki

Fortran cafe The Fortran Cafe'

Bethe-Salpeter wine The Bethe-Salpeter-Equation (BSE) wine

A street entitled to Yambo in Rome
A bar entitled to Yambo in Rome
A bar entitled to Yambo in Rome Yambo road, bar & restaurant

Basic strategy in computing surface spectra: optics and energy loss
by the Yambo Team

[Tutorial pdf document]

HREELS design by Lucia Caramella

In this tutorial you will learn how to carry out basic simulations of RAS and HREELS spectra at semiconducting surfaces using the Yambo code. In particular we will demonstrate how reasonable agreement with experimental data can be achieved even at the independent particle (RPA) level of approximation, and illustrate the effect of higher order corrections on the computed spectra. This feature is available only in the yambo_surf executable, so make sure you have compiled it!

The tutorial is split into two parts. In the first part we will deal with RAS (reflectance anisotropy spectroscopy) calculations on the GaAs(001)-c(4x4) surface using an asymmetric slab. In the second part we will consider a symmetric Si(001)-c(4x2) surface, demonstrating the calculation of RAS including many body effects, and the calculation of HREELS (high-resolution electron energy loss spectroscopy) data.

Tutorial structure

Once the tutorial zip file is unzipped the following folder structure will appear
COPYING README Surface_spectroscopy/
with the latter folder containing a directory for each structural model
>ls Surface_spectroscopy GaAs001-c4x4/ Si001-c4x2/
Each system in turn is composed of the subdirectories

Reflectance Anisotropy Spectroscopy: yambo_surf -s r

Theoretical background

The reflectance anisotropy signal is defined as the difference between normalized reflectivities for light polarized along perpendicular surface directions x and y. Since the reflectivity from the bulk (in cubic materials) is isotropic, this difference spectra is completely of surface origin. Equivalent expressions for the RAS were derived by Bagchi, Barrera and Rajagopal [5] and by McIntyre and Aspnes [3]. In the language of a three-layer model, the RAS is computed from:

Hence, the main quantity that is calculated in this runlevel is the dielectric tensor of the surface layer. We will mostly work within the independent particle scheme, in the long wavelength limit, and without local fields (i.e. equivalent to yambo -o c, with NGBlk=1). The tensor will be calculated for the (perpendicular) polarizations specified by q0x and q0y. Later we will show how to add some higher order effects.

The material: GaAs(001)-c(4x4)

  • FCC lattice
  • 7 layer slab, back surface (Ga) terminated by pseudo-hydrogens (Z=1.25)
  • Top layer composed of 3 buckled Ga-As dimers per cell.
  • 70 atoms per cell (236 electrons - 118 occupied bands)
  • Lattice constant 10.6 [a.u.]
  • Plane wave cutoff 13 Rydberg
GaAs structure

[01] Initialization: 01_init (yambo -i -V RL)

Let's begin with the GaAs(001)-c(4x4) surface. Enter the YAMBO folder
>cd GaAs001-c4x4/YAMBO/ >ls Experiment INPUTS REFERENCE SAVE
and run yambo_surf in the usual way to generate the setup databases:
yambo_surf -F INPUTS/01_init
Note that we use a fairly low number of G-vectors throughout just to keep things running quickly.

[01] Initialization: 01_init (yambo -i -V RL)

The RAS runlevel can be launched using yambo_surf -s r. Fairly complete descriptions of the input file variables are available on the RAS runlevel page. The pre-prepared input file INPUTS/02_RAS indicates the runlevel:
sursp # [R] Surface Spectroscopy ras # [R] Reflectance anisotropy spectroscopy (RAS)
Following some standard variables that you should be familiar with by now are variables specific to the RAS calculation.
% q0x 1.000000 | 0.000000 | 0.000000 | # [RAS] [cc] X-polarization direction % % q0y 0.000000 | 1.000000 | 0.000000 | # [RAS] [cc] Y-polarization direction %
These define the two surface polarization vectors in cartesian coordinates. They must be orthogonal, and form part of a right-handed frame with the surface normal.

There are also several lines allowing the bulk dielectric function to be read.
BulkFile= "Experiment/EPSILON_GaAs_bulk.dat" # [RAS] File containing bulk BulkForm= "2KK" # [RAS] Format of bulk (`3ri`,`3ir`,`2KK`,`3KK`) BlkBroad= 0.10000 eV # [RAS] Broadening of bulk eps
Have a quick look first at the contents of the bulk file, defined by BulkFile, and plot the first column with gnuplot. The form of the data in the bulk file is determined by BulkForm. Here we are using the experimental bulk epsilon of GaAs for no particular reason - one should really use a computed one as in the case of Si(001) below. As only the imaginary part is given, the real part is generated through a Kramers-Kronig transformation, using a broadening BlkBroad.

Basic output

Go ahead and launch the code using the predefined input file:
yambo_surf -F INPUTS/02_RAS_full ... <---> [04] Surface spectroscopy module <---> [04.01] Surface geometry and spectral analysis <---> [04.02] Bulk epsilon data import <---> Processed 101 data [ 0.00-10.00eV] from 42 records [ 2.03- 5.62eV] <---> [04.03] Calculating eps cell. <---> Absorption @ q || x : 0.70711 0.70711 0.00000 <---> [X-CG] R(p) Tot o/o(of R) : 18592 38704 100 <01s> Xo@q[1] 1-101 |####################| [100%] --(E) --(X) <01s> Absorption @ q || y : -0.70711 0.70711 0.00000 <02s> Xo@q[1] 1-101 |####################| [100%] --(E) --(X) <02s> [04.04] Write eps cell for plotting <02s> [04.05] Slab symmetries <02s> [04.06] Surface dielectric function <02s> [04.07] RAS <02s> [05] Game Over & Game summary
If you set up the bulk file incorrectly, the code will issue a warning and continue, thus computing anyway the SAVE/ndb.dipoles database. It's best not to use the -J option here, as the ndb.dipoles will be used again.

You should now have some new output files o.ras*, o.eps-cell*, o.eps-surf*. Check first inside the report file r_sursp_ras* that the geometry is as expected,
[04.01] Surface geometry and spectral analysis ============================================== Polarization directions: - X 1.000000 0.000000 0.000000 - Y 0.000000 1.000000 0.000000 Surface normal : 0.000000 0.000000 1.000000 Surface vector indices (ix,iy,iz) : 1 2 3 Supercell height az (a.u.) : 42.40000 No cut off function in matrix elements.
and now plot the contents of the RAS output file o.ras-rpa. Not so useful. Look at the diagram of the surface structure model. Do you see why?

Now change the polarization directions to correspond to meaningful directions in the surface plane: put q0x perpendicular to the dimer axes, q0z parallel to the dimer axes (the X and Y in the figure are the Cartesian axes, as used in the DFT input file) and run the code again.

The new output should look something like this:
[04.01] Surface geometry and spectral analysis ============================================== Polarization directions: - X 0.707107 0.707107 0.000000 - Y-0.707107 0.707107 0.000000
and the output file o.ras-rpa_01:
# Version 3.3.1 Revision 2087 # http://www.yambo-code.org # # 1 2 3 4 5 6 7 8 #hw (eV) RAS e1(x-y) e2(x-y) A B A x e2 B x e1 0.000000 0.000000 1.341996 0.000000 0.123842 0.000000 0.000000 0.000000 0.1000E-1 0.5830E-7 1.342 0.2266E-2 0.1238 0.1810E-4 0.2807E-3 0.2429E-4 0.2000E-1 0.2333E-6 1.342 0.4534E-2 0.1238 0.3613E-4 0.5615E-3 0.4850E-4
The different columns help in the interpretation of the computed RAS spectra, by splitting the RAS into surface and bulk components (note the bulk dielectric function in the denominator). Therefore, it is useful to rewrite the RAS as

This expression for the RAS makes it easier to identify structures in the RAS as being due to absorption (imaginary part) or dissipation (real part) in the surface layer [1]. The output of the yambo surface spectroscopy runlevel (o.ras-rpa) thus corresponds to the following:


You should now be able to identify which peaks in the computed RAS derive from the absorption, and thus can be easily analysed.

Real space cutoff

What we really want however, is the RAS of the front surface only, and not the response of the whole slab, including the passivated back surface. We would also like to relate spectral features to specific surface structures. To extract the front surface we use a real-space cutoff on the transition matrix elements:

The "zero" position and extent of the cutoff function are determined by the variables CutZero and CutStep . To use this, edit the relevant lines in the input, e.g. INPUTS/04_RAS_front (uncomment Cutoff):
% Cutoff # [RAS] Cutoff mode (Uncomment to use) CutZero= 0.250000 # [RAS] Zero position of cutoff fn (Frac) CutStep= 0.500000 # [RAS] Width of cutoff function (Frac)
Note that the values here are in fractional units of the cell axis, i.e. [0..1] range. Run the code now with
yambo_ras -F INPUTS/04_RAS_front -J 04_RAS_front
and note the creation of a new database in 04_RAS_front/ndb.dipoles_cut. This run is relatively slow because of a double sum on G. This will remove the contribution from the back surface. Let's go one more step and isolate the contribution from the topmost layers containing the dimers:
yambo_ras -F INPUTS/05_RAS_dimer -J 05_RAS_dimer
and have a look at the results all together.

We can understand three important things:

Last, let's see how our calculated spectrum compared with experiment. GaAs(001)-c(4x4) can be prepared in different ways. During MBE (molecular beam epitaxy) growth, a different reconstruction appears depending on whether As2 molecules (beta phase) or As4 molecules (alpha phase) are used. The measured RAS data appears in the Experiment folder: make a comparison with our "best" calculated spectrum. Can you now identify which phase corresponds to the Ga-As heterodimer structure we have assumed?

You're right - its the alpha phase. The beta phase in fact has pure unbuckled As-As dimers. Hence the surface richer in As (beta) is grown using As2 molecules, opposite to what you might expect!
The agreement is already very good, in spite of the approximations we have used in the calculation. What do you think is the most important parameter to converge/change in order to improve the agreement with experiment across the whole spectral range?


With the ypp_surf postprocessor one can analyse the peaks further by making a decomposition into different contributions (s-s, s-b, etc). In practice this is an essential part of the analysis, but is quite boring, and not worth studying further here.


  1. Use another value for the cutoff to extract a signal for the back layers. This should convince you that, even though the slab is electronically passivated, an optical anisotropy remains. (Hint: load the GaAs_layers.gnu gnuplot file in the PWSCF directory). Can you understand where this anisotropy comes from?

The material: Si(001)-c(4x2)

(2) Si(001)-c(4x2)
  • FCC lattice
  • 12 layer symmetric slab
  • Top/bottom layer composed of alternating buckled Si dimers.
  • 48 atoms per cell (192 electrons - 96 occupied bands)
  • Lattice constant 10.191 [a.u.]
  • Plane wave cutoff 13 Rydberg
Si structure

RAS for a symmetric slab

Let's move onto the Si(001)-c(4x2) surface. Enter the YAMBO folder
>cd ../Si001-c4x2/YAMBO/ >ls Experiment INPUTS REFERENCE SAVE
and again run yambo_surf to generate the setup databases:
yambo_surf -F INPUTS/01_init
Here we have a centrosymmetric slab, and don't necessarily need to use the cutoff function (do you see why a symmetric slab is not suitable for studying RAS of GaAs(001)?) You can quickly confirm this using:
yambo_surf -F INPUTS/02_RAS_full -J 02_RAS_full yambo_surf -F INPUTS/03_RAS_half -J 03_RAS_half
Note that here we are using a properly calculated Si bulk epsilon datafile (DFT-RPA), and have rigidly shifted it also by 0.5eV.
BlkShift= 0.500000 eV # [RAS] GW shift of bulk eps BulkFile= "Experiment/EPSILON_Si_bulk_RPA.dat" # [RAS] File containing bulk
Plot the two output files, and compare with experimental data in the Experiment folder.

Again, it already looks pretty good at the RPA level. This ab initio modeling is a bit of a joke!

Inclusion of higher order effects

A centrosymmetric slab allows us to move up a level, and add in some many body effects (at least until a real space cutoff is included...hi, Bernardo!). We do this through separate yambo calculations for the same two polarizations as before, and post process the results. Let's assume you've already done a GW calculation for the slab, which is no mean feat, and obtained an opening of the gap of +0.6eV. From this we can do a rough BSE calculation, just looking at the low energy part of the spectrum to save time. Remember the steps from previous lessons: first calculate the screening yambo -b, then construct and solve BSE kernel yambo -o b -k sex -y h.
yambo -F INPUTS/03_W yambo -F INPUTS/04_BSE_X yambo -F INPUTS/04_BSE_Y
Have a look inside these input files: you should be familiar with most of the variables by now. Compare the components of epsilon for the cell, computed within BSE and within RPA, from our earlier calculation.

Finally a trivial postprocessing step ypp_surf -a ras, where the two output files are read and the RAS computed in the normal way:
ypp_surf -F INPUTS/05_BSE_RAS
Just take care to define the parameter d correctly: it should be half the cell height in this case.

The result? Not a huge amount of difference between the RPA and the BSE results...It depends mostly on the value of the scissors operator correction. Can you understand why this is, by looking at the formula for the RAS?

Exercises (pick one or two!)

  • Use the cutoff function to extract a response for the fixed bulk layers of the slab. How does it compare with the bulk response calculated from the usual 2-atom cell?
  • Repeat the BSE calculation, but now using the Coulomb cutoff and random integration method (RIM) appropriate for a slab (2D sheet).
  • Reflectance EELS: yambo_surf -s e

    Theoretical basis

    Last we move on to the calculation of High Resolution EELS. In this experiment we measure the differential scattering cross section of electrons impinging on a surface, which can be written as:

    This expression consists of two parts: the kinematic term, which contains all the information about the experimental geometry, and the loss function, which contains all the interesting juicy physics of electron scattering. Within the dipole regime, we can use a true three-layer model to derive the loss function, similar to what was done for RAS. The most flexible and precise version of this model is that of Del Sole/Selloni[6]:


    So you see, also in the case of HREELS, the problem reduces to one of computing the dielectric function of a thin surface layer. The dependence on transferred momentum q is partially included, as we take the optical dielectric functions for convenience. This works quite well for typical experiments. All the tricks we discussed before regarding the cutoff function and so on apply also here, but be careful: the d parameter is actually used here, and must be determined correctly.

    Basic calculation

    The modeling is a little complicated due to the kinematic term and different possible scattering geometries. Here we just take you through the basic steps of the calculation.

    A typical input file is generated using yambo_surf -s e (extra options available with -V qp or -V gen):
    yambo_surf -s e -F INPUTS/06_REELS_full
    Inside you can see various new variables relating to the experimental geometry - see the various links from the REELS runlevel page for a more detailed explanation. Here we have to specific the thickness of the surface layer through the cutoff or the Layers variable. We've selected the loss function described above and a simple averaging over the detector window.

    Go ahead and run the code
    yambo_surf -F INPUTS/06_REELS_full -J 06_REELS_full
    and plot the output. You should be able to guess what the peaks at 1.8, 3.5, 5.1eV correspond to! (Hint: remember the RAS spectra).

    One major approximation being made here is that the electron does not penetrate the crystal surface. We can estimate its penetration depth from the universal scattering curve:

    At 40eV incident energy, it is about 6A, or 12 au. We can use this information then to augment the surface loss function with a term that represents excitation of bulk modes by the electron inside the material. This is done by defining a non-zero penetration depth:
    yambo_surf -s e -F INPUTS/06_REELS_bulk -J 06_REELS_bulk
    Run the code again, with this modified input file,
    yambo_surf -F INPUTS/06_REELS_bulk -J 06_REELS_bulk
    and now compare with the experimental data in the Experiment folder.

    It should be clear now that the higher energy features are excitations from the E1 and E2 critical points of the bulk.

    Exercises (pick one or two!)


    1. C. Hogan and R. Del Sole, Optical properties of the GaAs(001)-c(4x4) surface: direct analysis of the surface dielectric function, phys. stat. sol.(b) 242, 3040 (2005)
    2. C. Hogan, R. Del Sole, and G. Onida, Optical properties of real surfaces from microscopic calculations of the dielectric function of finite atomic slabs Phys. Rev. B 68 035405 (2003)
    3. J. D. E. McIntyre and D. E. Aspnes, Differential reflection spectro-scopy of very thin surface filmsSurf. Sci. 24 417 (1971)
    4. F. Manghi, R. Del Sole, A. Selloni, and E. Molinari, Anisotropy of surface optical properties from first-principles calculations Phys. Rev. B 41 9935 (1990).
    5. A. Bagchi, R. G. Barrera and A. K. Rajagopal, Perturbative approach to the calculation of the electric field near a metal surface, Phys. Rev. B 20 (1979)
    6. A. Selloni and R. Del Sole, Optical and electron energy-loss spectra of Si(111)-2x1 Surf Sci 168 35 (1986)
    7. R. Esquivel-Sirvent and Cecilia Noguez, Electron energy loss for anisotropic systems: Application to GaN (101bar0)Phys. Rev. B 587367 (1998)
    8. H. Luth, Surfaces and Interfaces of Solids Springer-Verlag, Berlin (1993)
    9. H. Ibach and D. L. Mills, Electron Energy Loss Spectroscopy and Surface Vibrations Academic Press, New York, (1982).