[Shortcut to examples 01 02 03 04 05 06 07 and exercises ]

by Yambo Team

*SiH4 structure.
*

- SiH
_{4}: common name Silane - Supercell 35.0 [a.u.]
- Plane waves cutoff 25 Ry (equal to 90000 RL vectors)

Silane - SiH_{4} 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.

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.

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 -Dand 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 RLand setting the MaxGvecs variable to a lower value.

Here we reduce the g-vectors to 50000, we initialize typing:

localhost:>yambo -F Inputs/01_init -J 01_initand 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

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:

InNGsBlkXd=300 RL # (Xd) Response block sizeand

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 300^{st}
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).

localhost:>yambo -F Inputs/02_alda_gspace -J 02_ALDA_g_spaceThe optional

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 fxc^{HEG} (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.

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_rspaceyou 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 summarythat 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 SiH_{4} molecule turns out to be a very challenging molecule as described in [2], where several functionals are tested displaying quite dissimilar results.

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__.

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:

% 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_couplingHere 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 BZgrows 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.

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 SiH_{4} 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.

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 dThen open the Inputs/07_BSE_diago file and uncomment the

WRbsWF # [BSS] Write to disk excitonic the FWsoption. 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

localhost:>ypp -e sand then put the hole in one of on one of the valence bands that more contribute to the electron-hole wavefunction.

- 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?

- SiH
_{4}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)