SiH4

WARNIGN: This page is presently under construction

The material

SiH4: common name Silane
Supercell 35.0 [a.u.]
Plane waves cutoff 25 Ry (equal to 90000 RL vectors)

Silane - SiH4 is a colorless gas with a repulsive odor at room temperature and atmospheric pressure. It is stable indefinitely in metal containers but in the pure state is spontaneously combustible in air.

Introduction

In this tutorial we will see how to calculate the optical spectrum of the Silane molecule, by means of the TDDFT, in the ALDA approximation, and the Bethe-Salpeter approach. The most common approach when treating localized systems, as atoms, cluster and molecules with plane-wave basis set, is the so-called supercell approach that consists in to put the system in a such large box that it can be considered isolated. At difference of bulk systems, now the dependence of the results on the supercell size, however, has to be checked carefully. Another big difference with respect bulk calculations is that no k-point sampling is needed for the wave function, so we can proceed considering only the gamma point.

In the following we will describe two different way to perform TDDFT/ALDA calculations. For each example you can find a link to the Input file with a brief explanation of the variables meaning.

We will also use this system to explain in detail some of the most relevant technical aspects of yambo.

Initialization

01_init (yambo -i -V RL)

This section describes the standard initialization procedure. As this initialization is common to all yambo runs it will be described extensively only in this case.

To run this example enter the YAMBO directory (where you have previously created the SAVE folder). Before running the initialization procedure let's have a look into the database produced by the yambo interfaces. In order to do that type:

localhost:>yambo -D 

and you will see

[RD./SAVE//ns.db1]------------------------------------------
 Bands                           :  40
 K-points                        : 1
 G-vectors             [RL space]:  90447
 Components       [wavefunctions]:  90447
 Symmetries       [spatial+T-rev]: 2
 Spinor components               : 1
 Spin polarizations              : 1
 Temperature                 [ev]: 0.000000
 Electrons                       : 8.000000
 WF G-vectors                    :  90447
 Max atoms/species               :  4
 No. of atom species             : 2

This is the header of the electronic structure database. It is interesting to note here the huge number of plane waves we have to work with. With respect a bulk calculation, now we have only the gamma point, but a big number of plane wave is needed to represent an isolated system. Of course this number increases with the supercell size. Usually it is not necessary to work with such a big number of plane waves that were needed for converging the ground state calculation. We can reduce this number typing:

yambo -i -V RL

and setting the MaxGvecs variable to a lower value. Note that this operation will limit the maximum number of wave-vector than can be used in all the calculations, and such approximation have to be carefully tested. Here we reduce the g-vectors to 50000, we initialize typing:

localhost:>yambo -F Inputs/01_init -J 01_init

and you will see

[02] CORE Variables Setup
 <---> [02.01] Unit cells
 <---> [02.02] Symmetries
 <---> [02.03] RL shells
 <---> Shells finder |                    | [000%] --(E) --(X)
 <---> Shells finder |####################| [100%] 01m-14s(E) 01m-14s(X)
 <---> [02.04] K-grid lattice
 <---> [02.05] Energies [ev] & Occupations
 <---> [03] Transferred momenta grid
 <---> [02.04] K-grid lattice
 <---> [02.05] Energies [ev] & Occupations
 <---> [03] Transferred momenta grid
 <---> X indexes |####################| [100%] --(E) --(X)
 <---> SE indexes |####################| [100%] --(E) --(X)
 <---> [04] Game Over & Game summary

If now you edit the r-01_setup you will find a lot of useless information about your system like a short summary of the electronic structure previously calculated at DFT level, and some details about atomic structure and the cell, for instance the HOMO-LUMO gap

 States summary         : Full        Metallic    Empty
                          0001-0004               0005-0040
 Indirect Gaps      [ev]: 7.914498  7.914498
 Direct Gaps        [ev]: 7.914498  7.914498
 X BZ K-points : 1

Adiabatic Local Density approximation (ALDA) in reciprocal space

02_ALDA (yambo -o c -t alda)

Now we can start for our TDDFT calculation in ALDA approximation. Here we do perform this calculation solving the TDDFT Dyson equation in G space representation: In Inputs/02_alda_gspace such a calculation is done with:

NGsBlkXd=300	RL	#	(Xd)	 Response block size

and

FxcGRLc= 300 RL    # [TDDFT] XC-kernel RL size 

where the response function size is set to 300 RL. This means including Local Field Effects corresponding to charge oscillations expanded up to the 300st RL component, and similarly also the matrix size of the ALDA kernel has been set to 300 RL. The converged value for NGsBlkXd and FxcGRLc shall be determined by doing several calculations with different values and checking the effect on the final physical result (the absorption spectrum, in this case). To run this example type

localhost:>yambo -F Inputs/02_alda_gspace -J 02_ALDA_g_space

The optional -J flag is used to label the output/report/log files (see here for the command line synopsis).

Then you can edit the input file changing the values of NGsBlkXd and FxcGRLc variables, for instance to 400 or more RL. Now we can compare the result of our calculations. The experimental result [1] denotes a rather broad spectrum where we can recognize three peaks located at: 8.8 eV, 9.7eV, 10.7eV. If we plot the obtained outputs:

The comparison with the experiment clearly shows that the agreement with the experiment is really poor. Most importantly, it shows that the absorption intensity is really small (if you enlarge the energy range you would observe that the spectrum here is hundreds of time too small). Such instability is due to a numerical instability of the TDDFT Dyson equation, with ALDA kernel, when applied to low-dimensional electronic systems, and evaluated in reciprocal space. The instability is due to the presence of large regions with vanishing density, as in large super-cells. In a region of the pace with vanishing density we have that fxcHEG (n) diverges so that the evaluation of the term involving the fxc in the Dyson equation cannot be directly calculated in reciprocal space, because the f ALDA is not well defined.

To overcome this problem yambo can solve the TDDFT equation within the ALDA in the basis of the e–h pairs, instead of the reciprocal space.

ALDA in transition space

03_ALDA_rspace(yambo -o b -k alda -y h)

Switching in the basis of the e-h pairs, when we have a static fxc (like in the ALDA) the TDDFT equation can be easily solved using the standard Bethe-Salpter techniques described in BSS. Now, the fxc kernel is evaluated in real space, and the problem reduces to the diagonalization of the Casida equations. Looking at the input file 03_ALDA_rspace we see that:

BSresKmod= "x"               # [BSK] Resonant Kernel mode. (`x`;`c`;`d`)
BScplKmod= "none"            # [BSK] Coupling Kernel mode. (`x`;`c`;`d`;`u`)

we are including only the resonant part (Tamm-Dancoff approximation).

BSENGexx=  50301       RL    # [BSK] Exchange components
BSSmod= "h"                  # [BSS] Solvers `h/d/i/t`

moreover we are including all the components of the local-fields and fxc kernel as it is treated in real-space. Furthermore we have chosen to use the Haydock recursive mode to diagonalize the matrix, which is particularly efficient with large matrices.

Running

localhost:>yambo -F Inputs/03_alda_rspace -J 03_alda_rspace

you will see

[05.01] Haydock solver
<--> [Haydock] Iteration 1
<--> [Haydock] Iteration 2
<--> [Haydock] Iteration 3 Accuracy :  0.02941| -0.02000
<--> [Haydock] Iteration 4 Accuracy :  0.07440| -0.02000
<--> [Haydock] Iteration 5 Accuracy :  0.02655| -0.02000
<--> [Haydock] Iteration 6 Accuracy :  0.05848| -0.02000
<--> [Haydock] Iteration 7 Accuracy :   0.0099| -0.02000
[06] Game Over & Game summary

that the procedure stops when the desired accuracy is reached.

Next let's remove the Tamm-Dancoff approximation, considering also the coupling between the resonant and anti-resonant part of the TDDFT matrix considering the input file 03_ALDA_rspace_cpl which differs from the previous one, for the flag:

BScplKmod= "x"            # [BSK] Coupling Kernel mode. (`x`;`c`;`d`;`u`)

that activates the coupling part to the matrix.

Plotting the output we get:

First, we note that, the inclusion of the coupling part in the matrix has a rather strong effects. When the Tamm-Dancoff approximation is removed we get the correct number of peaks in the energy range of interest, but unfortunately, even if now we get a stable result, it looks that ALDA does not reproduce the position of the experimental peaks. Note that such a failure is not general for finite system, as the ALDA is commonly used for computation of optical spectra of molecules and clusters [2]. The SiH4 molecule turns out to be a very challenging molecule as described in [2], where several functionals are tested displaying quite dissimilar results.

The Statically screened Electron-electron interaction

04_W(yambo -b)

As simple RPA and ALDA approximations have not described correctly the experimental spectrum we have to move towards more elaborate techniques like the Bethe-Salpeter (BS) equation.

A key ingredient in the BS kernel is the electron-electron interaction commonly evaluated in the static approximation. The input file Inputs/04_W describes how to calculate it. The variables used in this input file have the same physical meaning of those used in the optical absorption calculation. The only difference is that, in general, the response function dimension obtained in the examples (02-03), gives an upper bound to the number of RL vectors needed here. So let's type (no -J option here) ...

localhost:>yambo -F Inputs/04_W  -J 04_W

... and after waiting for some seconds we get the result, that is the database 04_W/ndb.em1s.

Bethe-Salpeter Equation (BSE)

05_bse(yambo -o b -k sex -y h)

The input file Inputs/05_BSE describes how to calculate an excitonic absorption spectrum, using as a solver of the BS equation the recursive Haydock Method. Before running yambo we must have a closer look at the input file: The line

% KfnQP_E
 5.10000 | 1.000000 | 1.000000 |      # [EXTQP BSK BSS] E parameters (c/v)
%

applies a 5.1 eV scissor in order to open the LDA gap.

BSEmod= "causal"             # [BSE] resonant/causal/x_coupling/xc_coupling

Here we are doing a BSE calculation including both exchange and correlation terms in the the resonant part, while the coupling part is set to zero. Still the resonant and anti-resonant parts are mixed using a causal receipe. This makes the BSE Hermitian, therefore the standard Haydock approach can be used.

% BSEBands
 1 |  38 |                   # [BSK] Bands range
%
BSENGBlk=  400            RL  # [BSK] Screened interaction block size
BSENGexx=  20000          RL  # [BSK] Exchange components

Here we specify the range of bands to be used and the RL components for the screened interaction (BSENGBlk=400) (read from the ndb.em1s file) and for the Hartree term (BSENGexx=20000). Notice that because we are interested only the absorption, the number of bands entering in the BSE can be smaller than the one needed to calculate the screened interaction W. This fact is important because the BS kernel is written in Bloch space and its size

BS kernel size = Valence Bands × Conduction Bands  × K-points in the whole BZ

grows very fast with system size or k-points number. In our case the size is 136. The BS bands range must be converged with care trying to keep as a few bands as possible.

Running

localhost:>yambo -F Inputs/05_BSE -J "05_BSE,04_W" &> output_yambo

will calculate the BS kernel and store in the 05_BSE folder the databases. Note the use of a double field for the -J option. In this way Yambo will seek for databases both in the 05_BSE and 04_W folders. Before analyzing the BS absorption spectrum we can edit the log file l-05_BSE_optics_bse_bss to observe the on-the-fly live timing feature of yambo:

[...]
<---> Kernel |                    | [000%] --(E) --(X)
<05s> Kernel |#                   | [007%] 05s(E) 01m-07s(X)
<11s> Kernel |###                 | [015%] 10s(E) 01m-06s(X)
[...]
<01m-02s> Kernel |##################  | [093%] 01m-01s(E) 01m-06s(X)
<01m-07s> Kernel |####################| [100%] 01m-06s(E) 01m-06s(X)
[...]

The two numbers at the end of each line represents the elapsed time(E) and the expected time(X). The latter estimates the time needed to complete that particular section (the BS kernel construction in this case).

Finally we can compare the resulting BS optical absorption with the experimental results and the independent particle(IP) approximation. As you can see the excitonic effects strongly modify the spectra respect to the IP, but the peaks are still not enough close to the experimental ones.

Now we will try to improve the outcome of our calculations including coupling terms in the BSE.

BSE with coupling terms

06_bse (yambo -o b -k sex -y h)

We remove the Tamm-Dancoff approximation, as we have done in the ALDA case including the coupling between the resonant and anti-resonant part of the BSE matrix. The new input file 06_BSE_cpl which differs from the previous one just for the flag:

BSEmod= "xc_coupling"        # [BSE] resonant/causal/x_coupling/xc_coupling

that turns on the coupling part of the BSE. Running the calculation:

localhost:>yambo -F Inputs/06_BSE_cpl -J "06_BSE_cpl,04_W"

you will immediately notice that the size of the BS Kernel increased from 136 to 272 due to the coupling terms. Moreover the full BS matrix is now not Hermitian, but also in this case it is possible to use the Haydock technique thanks to a recent generalization to pseudo-Hermitian matrices[6]. Now we compare the result with the previous one:

As you can see including coupling terms improve the agreement with the first ant third peak, but the second one remains unchanged. Why BSE is not able to reproduce the experimental result? This failure stems from the fact that we used a rigid scissor operator for all the bands, and this approximation is usually not verified in molecular systems where the GW corrections can depend from the character and the localization of a given orbitals. Moreover it has been shown in reference [4] that the usual diagonal approximation for the self-energy in the calculation of the quasi-particle correction fails for the SiH4 molecule. We advise to read references [3,4,5] in order to have an idea of the different methods and approximations employed to calculate the spectra of this "simple" molecule.

Analyzing the results with ypp

In the YAMBO package is included a very useful tool ypp to postprocess the different calculations and analyze the results. In this section we will use tool to analyze the peaks obtained in the previous calculations. In order to do so, we need to fully diagonalize the BS Hamiltonian, therefore repeat the previous calculations but choose as solver "d", diagonalization.

localhost:>yambo_gpl -F Inputs/06_BSE -o b -k sex -y d

Then open the Inputs/07_BSE_diago file and uncomment the

WRbsWF                       # [BSS] Write to disk excitonic the FWs

option. Run the code.

localhost:>yambo -F Inputs/06_BSE -J "07_BSE_diago,04_W"

Then we can order the peaks in the absorption according to their intensity and energy, in such a way to select the peak we are interested in, running:

localhost:>ypp -e s -J 07_BSE_diago

ypp will produce two files o-07_BSE_diago.exc_I_sorted and o-07_BSE_diago.exc_E_sorted . In the first the peaks of the BSE are ordered by intensity and in the second by energy. Now we want to analize the first one at E=8.6 eV, the number 111. ypp allows you to plot the excitonic wave-function. The electron-hole is a two particle wave-function and tells you where it is possible to find the electron for a give position of the hole. Now the point is, where to put the hole? The best way to put to put the hole is to decompose the excitonic peak in the basis of free electron-hole excitations with the comand:

localhost:>ypp -e s

and then put the hole in one of on one of the valence bands that more contribute to the electron-hole wavefunction.

Additional Exercises

• Plot the spectrum obtained in section 03 and check the effect of increasing the polarization bands.
• Calculate the absorption spectra in the RPA approximation, including Local Fields, both in reciprocal space and in the basis of e-h pairs: 02-03 and check the effect of increasing the block size of the response function and evaluate the effect of the ALDA kernel with respect the local field effects.
• Perform again BSE calculations 05 but this time include only exchange terms or diagonal terms only in the Coupling and Resonant part of the BSE, how much the result changes?

References

• SiH4 on wikipedia.
• MAL Marques, A. Castro and A Rubio, J. Chem. Phys. 115,3006 (2001).
• U. Itoh, Y. Toyoshima, H. Onuki, N. Washida, and T. Ibuki3, J. Chem. Phys. 85, 4867 (1986)
• J. C. Grossman, M. Rohlfing, L. Mitas, S. G. Louie, and M. L. Cohen, Phys. Rev. Lett. 86, 472–475 (2001)
• R. J. Needs, P. R. C. Kent, A. R. Porter,. M. D. Towler, G. Rajagopal, Int. Journal of Quantum Chem.,86, 218 (2002)
• M. Grüning, A. Marini and X. Gonze, Nano Letters,9, 2820 (2009)