Hydrogen chain

= Fun with a Hydrogen chain:
the obscure(?) reasons for the failure of the ALDA
by Andrea Marini =

[Tutorial [[../files/hydrogen_chain.pdf|pdf]] document]

In the first TDDFT [[../tddft/index.php|tutorial]] we learned how to calculate the response function in the TDDFT scheme.

You should have noticed that the performance of the ALDA gradually worsens moving from 0 to 3 dimensions. This means that the drawbacks of the approximation are somehow due to possibility that the electrons move in "wider" regions of space.
The subject of this tutorial is to show how yambo can be used to pin down the reasons for the failure of the ALDA in a simple system. The conclusions of this tutorial do not mean to be general, but they should convince you that there is a link between the local assumption of the ALDA and the polarization of electrons in extended directions.

The system used in this tutorial is an infinite H₂ molecular chain. This is a very simple physical system consisting of H atoms that are distributed in sets of two atoms placed at a variable distance X from each other.

The physical properties of the chain are functions of the distance X. When X=2. a.u. the system is metallic. By increasing X the chain becomes semiconducting with increasing gap.

First of all, after having downloaded the zip files of the tutorial yambo databases you should have six folders corresponding to the six values of X (2.05, 2.1, 2.2., 2.3, 2.4 and 2.5).

The ALDA failure

The first step of this tutorial is to run the calculation of the dynamical absorption in the TDLDA approximation.
We will describe here only the case of X=2.05 a.u.. You are invited to repeat the calculation for X=2.5 and, if you wish, for all the other cases.

Enter the 2.05 directory and run the setup launching yambo

localhost> cd 2.05
localhost>ls
BSE/  SAVE/
localhost> yambo

Editing the r_setup file we notice that the system has a small (0.27 eV) gap.
Now we can directly run the TDLDA calculation by typing

localhost>  yambo -o b -t a -y d -V qp

You are now redirected to the editing of the yambo.in input file.

optics                       # [R OPT] Optics
bse                          # [R BSK] Bethe Salpeter Equation.
alda_fxc                     # [R TDDFT] The ALDA TDDFT kernel
bss                          # [R BSS] Bethe Salpeter Equation solver
KfnQPdb= "none"              # [EXTQP BSK BSS] Database
KfnQP_N= 1                   # [EXTQP BSK BSS] Interpolation neighbours
% KfnQP_E
 0.000000 | 1.000000 | 1.000000 |      # [EXTQP BSK BSS] E parameters (c/v)
%
% KfnQP_W
 0.000    | 0.000    | 0.000    | 0.000    |     # [EXTQP BSK BSS] W parameters  (c/v)
%
KfnQP_Z= ( 1.000000 , 0.000000 )       # [EXTQP BSK BSS] Z factor  (c/v)
BSresKmod= "x"               # [BSK] Resonant Kernel mode. (`x`;`c`;`d`)
% BSEBands
  1 | 20 |                   # [BSK] Bands range
%
BSENGexx=  7659        RL    # [BSK] Exchange components
BSSmod= "d"                  # [BSS] Solvers `h/d/i/t`
% BEnRange
  0.00000 | 10.00000 | eV    # [BSS] Energy range
%
% BDmRange
  0.10000 |  0.10000 | eV    # [BSS] Damping range
%
BEnSteps= 100                # [BSS] Energy steps
% BLongDir
 1.000000 | 0.000000 | 0.000000 |      # [BSS] [cc] Electric Field
%

Please change the highlighted values to ...

% KfnQP_E
 3.500000 | 1.000000 | 1.000000 |      # [EXTQP BSK BSS] E parameters (c/v)
%
...
% BSEBands
  1 | 2 |
%
...
BEnSteps= 1000                # [BSS] Energy steps
...

Now use the resp verbosity to activate the flag needed to dump to file the eigenvectors of the BS Hamiltonian:

localhost>  yambo -o b -t a -y d -V resp

and remove the # to the flag

#WRbsWF                      # [BSS] Write to disk excitonic the FWs

Now run yambo. In the folder ./BSE you will find the absorption spectra (o-BSE.eps_q001-bd) calculated using the Many-Body based Bethe-Salpeter(BS) equation. If you wish you can find here a dedicated section about the BS calculation. For the moment we will keep the results of the Bethe-Salpeter equation as a reference for our calculations. The BS equation leads indeed a proper and accurate description of the optical properties of the H₂ molecular chains.
If you plot the ALDA result against the BS one you will find that there is a reasonable agreement between the two curves.

[[File:chain_2.05_alda_vs_bse.png ]]Now, please, repeat the same procedure in the 2.5 folder

localhost> cd ../2.5
localhost> yambo
localhost>  yambo -o b -t a -y d -V resp
...

If you plot again the ALDA result against the BS one you will find that in this case the performance of the ALDA is much worse.

[[File:chain_> ]]
The question now is ...

A closer look: plots

At this stage we need to get some more information from the TDLDA calculations that can be compared with the "exact" BS results. Let's proceed by plotting the electronic density and the wavefunction corresponding to the most intense peak in the dynamical polarizability.

A closer look: plots

The electronic density

The easiest quantity we can compare for the two chains (2.05 and 2.5 intratomic distance) is the electronic density.
To this end we can use the yambo post/pre processor (YPP). YPP is capable of performing some basic analyses using the pre-calculated informations stored in the yambo databases.
YPP works like yambo. It uses a series of options to create an input file containing only the variables that are relevant to that type of calculation.

localhost>cd 2.05/
localhost>ypp -H
 ___ __  _____  __ __  _____   _____ 
|   Y  ||  _  ||  Y  ||  _  \ |  _  |
|   |  ||. |  ||.    ||. |  / |. |  |
 \   _/ |. _  ||.\ / ||. _  \ |. |  |
  |: |  |: |  ||: |  ||: |   \|: |  |
  |::|  |:.|:.||:.|:.||::.   /|::.  |
  `--"  `-- --"`-- --"`-----" `-----"

 Tool: ypp 3.2.1 rev.506
 Description: Y(ambo) P(ost) P(rocessor) 

 -h :Short Help
 -H :Long Help
 -J    :Job string identifier
 -V    :Input file verbosity [opt=gen]
 -F    :Input file
 -I    :Core I/O directory
 -O    :Additional I/O directory
 -C    :Communications I/O directory
 -N :Skip MPI initialization
 -S :DataBases fragmentation
 -k    :BZ Grid generator [(k)pt,(q)pt,(l)ongitudinal]
 -e    :Excitons  [(s)ort,(a)mplitude,(w)ave]
 -s    :Electrons [(w)ave,(d)ensity,do(s)]
 -f :Free hole position [excitons plot]
 -r :BZ energy RIM analyzer

By      yambo developers group
        http://www.yambo-code.org

To plot the density we need to type

localhost> ypp -p d

The input file ypp.in will be

density                      # [R] Density
electrons                    # [R] Electrons (and holes)
Format= "g"                  # Output format [(g)nuplot/(x)crysden]
Direction= "1"               # [rlu] [1/2/3] for 1d or [12/13/23] for 2d [123] for 3D
FFTGvecs=  7659          RL  # [FFT] Plane-waves

As our system is one-dimensional and it is lying parallel to the z axis we need to change the Direction="12" so to perform a contour plot on the "XY" plane. To use Xcrysden switch to "x" the value of the Format variable.
Now you can launch ypp.
At the end of the quick calculation you will find in your directory the file o.density_2d.xsf that you can plot using Xcrysden. We have provided a script to automatically view this file in the ../bin folder. To use it type

localhost>../bin/launch_xcrysden.sh o.density_2d.xsf

Follow the same procedure for the other distance and compare the two densities. You should see something like

File:Chain 2.05.png

We see that the smaller the gap (@ 2.05 a.u.) the more delocalized the electronic density. The consequence is that electrons can move more freely in the 2.05 case. At the same time, however, the polarization will be more "metallic-like" where the ALDA is expected to work better.

The TDLDA excitations wavefunction

Using YPP we can do something more. We can plot the wavefunction corresponding to the most intense peak in the dielectric function. The excitation wavefunction is a two-coordinate (r₁, and r₂) function and represents the probability amplitude of finding the electron at position r₁ when the hole is in r₂. The larger this probability the more delocalized is the perturbation induced by the external field.
A mean large distance between the electron and hole also reflects the poor correlation felt by the two bodies. YPP can plot either the case where the electron and the hole are forced to be in the same space point, or the case where the hole's position is fixed somewhere and the electron is left free to move. This second case is more appropriate here.
We need first to identify the index of the corresponding eigenstate.
Running

localhost>ypp -e s

YPP will create a file named o.exc_I_sorted that contains the list of peaks in the dielectric function ordered with increasing peak intensity.

#
#  E [ev]     Strength   Index
#
  3.783688   1.000000   1.000000
   4.38457    0.03179    3.00000
  5.171     0.3731E-2   5.000

We see that the peak with energy 3.78 dominates the spectra and it corresponds to the index 1.
If we type now

localhost>ypp -e w

to get the input file

exc_wf                       # [R] Excitonic Wavefunction
excp                         # [R] Excitonic Properties
plot                         # [R] Plot
Format= "x"                  # Output format [(g)nuplot/(x)crysden]
Direction= "12"             # [rlu] [1/2/3] for 1d or [12/13/23] for 2d [123] for 3D
FFTGvecs=  7659          RL  # [FFT] Plane-waves
States= "1 - 1"              # Index of the BS state(s)
Degen_Step=   0.0100     eV  # Maximum energy separation of two degenerate states
% Cells
 1 |  1 |  1 |                        # Number of cell repetitions (even or 1)
%
% Hole
 0.000    | 0.000    | 0.000    |      # [cc] Hole position in unit cell
%

It is interesting to do a 3D plot. To this end we set Direction="123" and replace 1 with 16 in Cells (this value will expand the unit cell along the X direction). The States variable is already set to the peak number 1.
We put the hole in the middle of the chain by setting

% Hole
 1.02500 | 12.50000 | 12.50000  |           # [cc] Hole position in unit ce
%

remember that the first field must be set to the interatomic distance divided by two. In the case of 2.5 a.u. the first field will be 1.25.
If we run now ypp it will create the file o.exc_3d_1.xsf that you can plot by using

localhost>../bin/launch_xcrysden.sh o.density_2d.xsf

If we compare the plot with the one obtained by solving the BS equation we should see

Excitation wavefunctions for 2.05 a.u. distance

[[File:chain_2.05.o-BSE_RIM.exc_3d_1.xsf.png ]][[File:chain_2.05.o-TDLDA.exc_3d_1.xsf.png ]]|-

BS equation

TDLDA

Excitation wavefunctions for 2.5 a.u. distance

[[File:chain_2.5.o-BSE_RIM.exc_3d_1.xsf.png ]][[File:chain_2.5.o-TDLDA.exc_3d_3.xsf.png ]]|-

BS equation

TDLDA

We immediately see that, while the ALDA excitation is always spread all over the chain, in the 2.5 distance case the "true" wavefunction acquires a tail that decreases the probability of finding the electron and the hole very far from each other.
This saturation of the excitation wavefunction is the physical reason for the poor performance of the ALDA. As it is a local approximation it cannot take into account a long-range correlation between the electron and the hole.

A long-range kernel beyond the TDLDA

So, we realized the a simple plot of the excitation wavefunction can pin down the possible reasons for the breakdown of the ALDA. In general if you know the problem, you should be half way through the quest for a solution. Is it true?

The Bethe-Salpeter equation: tricks and tips of the 1D systems

The tricky case of the H₂ chain. When the BS equation can be (even) more complicated!

Additional Exercises

Repeat the whole tutorial for the other intratomic distances, whose databases are provided in the tutorial zip file.