SiH4: isolated molecule
[Shortcut to examples
01
02
03
04
05
06
and exercises
]
by Yambo Team
[Tutorial pdf document]
SiH4 structure.
The material
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.
[01] 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 -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. 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_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
[02] Adiabatic Local Density approximation (ALDA) 02_ALDA (yambo -o c -t a)
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 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 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 typelocalhost:>yambo -F Inputs/02_alda_gspace -J 02_ALDA_g_spaceThe 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 results are rather unstable. 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.
[03] ALDA: 03_ALDA_rspace(yambo -o b -t a -y h -V resp)
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`We are considering all the components of the local-fields and fxc kernel as it is treated in real-space. We have chosen to use the Haydock recursive mode to diagonalize the matrix, which is particularly efficient for large matrices. The variable
BSHayTrs= 0.001000define the threshold of the Haydock iterative procedure.
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.280702|0.1000E-2
[Haydock] Iteration 4 Accuracy : 0.390927|0.1000E-2
...
[Haydock] Iteration 30 Accuracy : 0.0017|0.1000E-2
[Haydock] Iteration 31 Accuracy : 0.0011|0.1000E-2
[Haydock] Iteration 32 Accuracy :0.9733E-3|0.1000E-2
[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:
We note that the inclusion of the coupling part in the matrix modifies the position and intensity of the main peak close to 11 eV bringing it closer to the experimetal value. Nevertheless in general 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.
[04] The Statically screened Electron-electron interaction: 04_W(yambo -b)
As simple RPA and ALDA 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
... and after waiting for some seconds we get the result, that is the database SAVE/ndb.em1s. Here it is better to skip the -J option because we will need to read this database in the next examples.
[05] BSE: 05_Bethe-Salpeter Equation(yambo -o b -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.
BSresKmod= "xc" # [BSK] Resonant Kernel mode. (`x`;`c`;`d`) BScplKmod= "none" # [BSK] Coupling Kernel mode. (`x`;`c`;`d`;`u`)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. 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 &> output_yambo
will calculate the BS kernel and store in the 05_BSE folder the databases. 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:
[...] <---> BSK | | [000%] --(E) --(X) <05s> BSK |# | [007%] 05s(E) 01m-07s(X) <11s> BSK |### | [015%] 10s(E) 01m-06s(X) [...] <01m-02s> BSK |################## | [093%] 01m-01s(E) 01m-06s(X) <01m-07s> BSK |####################| [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 IPA (with scissor operator). The situation in improved with respect of the TD-LDA. Indeed in the region of interest BSE shows three peaks, though their position does not agree very well with the experimental one.
Now we will including the coupling terms in the BSE as we have done previously for ALDA.
The new input file 05_BSE_cpl differs from the previous one just for the flag:
BScplKmod= "xc" # [BSK] Coupling Kernel mode. (`x`;`c`;`d`)
that turns on the coupling part of the BSE. Running the calculation:
localhost:>yambo -F Inputs/05_BSE_cpl -J 05_BSE_cplyou will 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]. Including the coupling modifies the spectrum with respect to the Tamm-Dancoff:
As in the ALDA the effect of including the coupling is to red-shift and reduce the intensity of some of the peaks (this is general for confined systems). In particular this is the case for the second and third peak. The latter is now on top of the experimental result, still the position of the first/second peaks is not very well described.
Why BSE is not accurately reproducing the experimental result?
The disagreement stems most probably 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 may depend strongly on the character and the localization of a given orbitals.
Furthermore 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.
[06] Analyzing the results with ypp
In the YAMBO package it is included a very useful tool ypp to postprocess the different calculations and analyze the results. In this section we will use this tool to analyze the peaks obtained in the previous calculations.
So let's start again from the input 05_Bethe-Salpeter Equation(file 05_BSE) and analyze the result. In order to do so, we need to fully diagonalize the BS Hamiltonian, therefore we repeat the calculation but we choose as solver "d", equal to diagonalization, and then we add the flag WRbsWF (verbosity -V resp) to calculate the BS eigenvectors.
Now we can order the peaks in the absorption spectra according to their intensity and energy, in such a way to select the one we are interested in, running:
localhost:>ypp -e s -J 05_BSEypp will produce two files o.exc_I_sorted and o.exc_E_sorted . In the first the peaks are ordered by intensity while in the second by energy. Now we want to analyze the one at E=8.604 eV, number 6.
ypp allows you to plot the excitonic wave-function. This wave-function is a two particle function that 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 the hole is to decompose the excitonic peak in the basis of free electron-hole excitations with the command:
localhost:>ypp -e a -J 05_BSEypp produces two files, in the first o.exc_amplitude_at_6 there the amplitudes of the different electron-hole excitation ordered according to their energy, and in the second o.exc_weight_at_6 the most important contribution to the exciton. Now you can plot most relevant valence bands with the command ypp -s w. At this point you will have an idea where to put the hole for the final plot of the excitonic wave-function with the command ypp -e w -J 05_BSE.
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.
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)