Time-dependent density functional theory (standard kernel): Difference between revisions

From The Yambo Project
Jump to navigation Jump to search
 
(16 intermediate revisions by one other user not shown)
Line 1: Line 1:
In this tutorial you will learn how to calculate optical spectra at the time-dependent density-functional theory level for bulk hBN.
In this module you will learn how to calculate optical spectra at the time-dependent density-functional theory level for bulk hBN.
== Background ==
== Background ==
The macroscopic dielectric function is obtained by including the so-called local field effects (LFE) in the calculation of the response function. Within the time-dependent DFT formalism this is achieved by solving the Dyson equation for the susceptibility ''X''. In reciprocal space this is given by:
Within the time-dependent DFT formalism the macroscopic dielectric function [[Local_fields|is obtained]] by solving the Dyson equation for the susceptibility ''X''. In reciprocal space this is given by:
[[File:Yambo-CH7.png|none|x50px|Yambo tutorial image]]
[[File:Yambo-CH7.png|none|x50px|Yambo tutorial image]]
The microscopic dielectric function is related to ''X'' by:
[[File:Yambo-CH5.png|none|x30px|Yambo tutorial image]]
and the macroscopic dielectric function is obtained by taking the (0,0) component of the inverse microscopic one:
[[File:Yambo-CH6.png|none|x50px|Yambo tutorial image]]
Experimental observables like the optical absorption and the electron energy loss can be obtained from the macroscopic dielectric function:
[[File:Yambo-CH8.png|none|x50px|Yambo tutorial image]]


The ''f <sup> xc</sup>'' term is called the exchange-correlation (xc) kernel and describes the response of the xc potential at a time ''t'' to changes in the density at all previous times. Approximations   
The ''f <sup> xc</sup>'' term is called the exchange-correlation (xc) kernel and describes the response of the xc potential at a time ''t'' to changes in the density at all previous times. Approximations   
Line 15: Line 9:
Within this family we consider here the adiabatic local density approximation (ALDA)
Within this family we consider here the adiabatic local density approximation (ALDA)
[[File:Eq_ALDA.png|none|x50px]]
[[File:Eq_ALDA.png|none|x50px]]
As this module will show for bulk hBN this approximation provides spectra similar to the RPA when applied to extended systems (it has instead a significant impact on molecular excitations).
One of the main problems that has been pinpointed in the literature is that ALDA and its likes miss a long range contribution essential to reproduce excitations in solids. Several approximations have been proposed  to introduce this essential long range contribution.<ref>L. Reining, V. Olevano, A. Rubio, and G. Onida, Phys. Rev. Lett. 88, 066404 (2002)</ref>,<ref>P. E. Trevisanutto, A. Terentjevs, L. A. Constantin, V. Olevano, and F. Della Sala, Physical Review B 87, 205143 (2013)</ref><ref>J. A. Berger, Phys. Rev. Lett. 115, 137402 (2015)</ref><ref>S. Rigamonti, S. Botti, V. Veniard, C. Draxl, L. Reining, and F. Sottile, Phys. Rev. Lett. 114, 146402 (2015)</ref>,<ref>S. Sharma, J. K. Dewhurst, A. Sanna, and E. K. U. Gross, Phys. Rev. Lett. 107, 186401 (2011)</ref>
Within this family we consider here the long-range correction approximation (LRC):
[[File:Eq_LRC.png|none|x50px]]
where &alpha; is an empirical parameter related to the system screening.<ref>S. Botti, F. Sottile, N. Vast, V. Olevano, L. Reining, H.-C.Weissker, A. Rubio, G. Onida, R. Del Sole, and R. Godby, Physical Review B 69, 155112 (2004)</ref>
Note that approximations of this form cannot be rigorously justified within the sole time-dependent density-functional framework. Additional key quantities such as the electron density current or the macroscopic polarization need to be considered.<ref>M. Gr&uml;ning, D. Sangalli, C.  Attaccalite, C. Physical Review B '''94''', 035149 (2016)</ref>,<ref>P. L. de Boeij, F. Kootstra, J. A. Berger, R. van Leeuwen, and J. G. Snijders, The Journal of Chemical Physics 115, 1995 (2001)</ref><ref>J. A. Berger, Phys. Rev. Lett. 115, 137402 (2015)</ref>


== Prerequisites ==
== Prerequisites ==
 
[[File:Yambo-handbook-v4.1.2-p-5.png|thumb|Cheatsheet on TDDFT|150px]]
* You must first complete the "How to use Yambo" modules
* You must first complete the "How to use Yambo" and the "Optics at the independent particle level" and "local fields" modules  


'''You will need''':
'''You will need''':
Line 35: Line 20:


== Choosing input parameters ==
== Choosing input parameters ==
Enter the folder for bulk hBN that contains the ''SAVE'' directory, run the initialization and generate the input file. From <code>yambo -H</code> you should understand that the correct option is <code>yambo -o c</code>. Let's add some command line options:
Enter the folder for bulk hBN that contains the SAVE directory, and generate the input file. From <code>yambo -H</code> you should understand that the correct option is <code>yambo -o c -k <opt> </code>.
<!--If you have already worked in this folder, you might like to clean up any old files, and run the initialization again:
Where <opt> can be 'hartree' as in the tutorial on the [[Local fields|local fields]], or 'alda'  or 'lrc' ([[Time-dependent density functional theory (long-range corrected kernel)|see next tutorial]]) respectively for the ALDA and LRC approximations for the xc kernel outlined above.
  cd YAMBO_TUTORIALS/hBN/YAMBO
Let's start by generating the input for a calculation with the ALDA kernel:
  rm yambo.in SAVE/ndb.*
  $ cd YAMBO_TUTORIALS/hBN-3D/YAMBO
  yambo  
  $ yambo       ''(Initialization)''
Take care to '''not''' delete the <code>ns.*</code> files. -->
  $ yambo -F yambo.in_ALDA -J 3D_ALDA -o c -k alda
  $ cd YAMBO_TUTORIALS/hBN/YAMBO
We thus use a new input file ''yambo.in_ALDA'', switch on the <code>FFTGvecs</code> variable, and label all outputs/databases with a ''3D_ALDA'' tag. Make sure to set/modify all of the following variables to:
  $ yambo              ''(initialization)''
  [[Variables#Chimod|Chimod]]= "ALDA"            # [X] IP/Hartree/ALDA/LRC/BSfxc
  $ yambo -F yambo.in_IP -J Full -o c
  [[Variables#FxcGRLc|FxcGRLc]]= '''   3        Ry'''   # [Xd] Response block size
This corresponds to optical properties in G-space at the independent particle level: in the input file this is indicated by (<code>Chimod= "IP"</code>).
  [[Variables#NGsBlkXd|NGsBlkXd]]=    3        Ry    # [Xd] Response block size
% [[Variables#QpntsRXd|QpntsRXd]]
  1 |  1 |                  # [Xd] Transferred momenta
%
% [[Variables#EnRngeXd|EnRngeXd]]
  0.00000 | 10.00000 | eV    # [Xd] Energy range
%
% [[Variables#DmRngeXd|DmRngeXd]]
0.200000 | 0.200000 | eV    # [Xd] Damping range
%
[[Variables#ETStpsXd|ETStpsXd]]= 1001              # [Xd] Total Energy steps
% [[Variables#LongDrXd|LongDrXd]]
1.000000 | 1.000000 | 0.000000 |        # [Xd] [cc] Electric Field
%
In this input file, we have selected:
* the long-wavelength limit q = 0 (optical limit)
* A ''q'' parallel to the BN planes
* The ALDA kernel, and expanded G-vectors in the screening and in the xc kernel up to 3 Ry (about 85 G-vectors)


==Optics runlevel==
==Optics runlevel==
For optical properties we are interested just in the long-wavelength limit ''q = 0''. This always corresponds to the ''first'' ''q''-point in the set of possible ''q=k-k' ''-points.
Change the following variables in the input file to:
% [[Variables#QpntsRX|QpntsRXd]]
  1 |  '''1''' |                  # [Xd] Transferred momenta
%
[[Variables#ETStpsX|ETStpsXd]]= '''1001'''              # [Xd] Total Energy steps
in order to select just the first ''q''. The last variable ensures we generate a smooth spectrum.
Save the input file and launch the code, keeping the command line options as before (i.e., just remove the lower case options):
Save the input file and launch the code, keeping the command line options as before (i.e., just remove the lower case options):
  $ yambo -F yambo.in_IP -J Full
  $ yambo -F yambo.in_ALDA -J 3D_ALDA
  ...
  ...
  <---> [05] Optics
  <01s> [05] Optics
  <---> [LA] SERIAL linear algebra
  <01s> [LA] SERIAL linear algebra
  <---> [DIP] Checking dipoles header
  <01s> [WF] Performing Wave-Functions I/O from ./SAVE
<---> [x,Vnl] computed using 4 projectors
  <01s> [FFT-Rho] Mesh size: 9   9  21
<---> [M  0.017 Gb] Alloc WF ( 0.016)
  <01s> [xc] Functional Slater exchange(X)+Perdew & Zunger(C)
<---> [WF] Performing Wave-Functions I/O from ./SAVE
  <01s> [xc] LIBXC used to calculate xc functional
<01s> Dipoles: P and iR (T): |########################################| [100%] 01s(E) 01s(X)
<01s> [M  0.001 Gb] Free WF ( 0.016)
  <01s> [DIP] Writing dipoles header
<01s> [X-CG] R(p) Tot o/o(of R) :   5501  52992    100
  <01s> Xo@q[1] |########################################| [100%] --(E) --(X)
  <01s> [06] Game Over & Game summary
 
$ ls
Full  SAVE  yambo.in_IP  r_setup
o-Full.eel_q1_ip  o-Full.eps_q1_ip  r-Full_optics_chi
Let's take a moment to understand what Yambo has done inside the Optics runlevel:
# Compute the ''[x,Vnl]'' term
# Read the wavefunctions from disc [WF]
# Compute the ''dipoles'', i.e. matrix elements of '''p'''
# Write the dipoles to disk as ''Full/ndb.dip*'' databases. This you can see in the report file:
$ grep -A20 "WR" r-Full_optics_chi
[WR./Full//ndb.dip_iR_and_P]
Brillouin Zone Q/K grids (IBZ/BZ):  14  72  14  72
RL vectors                  (WF): 1491
Electronic Temperature        [K]: 0.0000000
Bosonic    Temperature        [K]: 0.0000000
X band range          :  1  100
RL vectors in the sum  : 1491
[r,Vnl] included      :yes
  ...
  ...
<ol start="5">
<li>Finally, Yambo computes the non-interacting susceptibility ''X<sub>0</sub>'' for this ''q'', and writes the dielectric function inside the ''o-Full.eps_q1_ip'' file for plotting</li>
</ol>
==Energy cut off==
Before plotting the output, let's change a few more variables. The previous calculation used ''all'' the G-vectors in expanding the wavefunctions, 1491. This corresponds roughly to the cut off energy of 40Ry we used in the DFT calculation. Generally, however, we can use a smaller value. We use the verbosity to switch on this variable, and a new ''-J'' flag to avoid reading the previous database:
$ yambo -F yambo.in_IP '''-J 6Ry -V RL''' -o c
Change the '''value''' of <code>[[Variables#FFTGvecs|FFTGvecs]]</code> and also its '''unit''' from <code>RL</code> (number of G-vectors) to <code>Ry</code> (energy in Rydberg):
[[Variables#FFTGvecs|FFTGvecs]]= '''6'''          '''Ry'''    # [FFT] Plane-waves
Save the input file and launch the code again:
  $ yambo -F yambo.in_IP '''-J 6Ry -V RL'''
and then plot the ''o-Full.eps_q1_ip'' and ''o-6Ry.eps_q1_ip'' files:
$ gnuplot
gnuplot> plot "o-Full.eps_q1_ip" w l,"o-6Ry.eps_q1_ip" w p


[[File:CH-hBN-6Ry.png|none|500px|Yambo tutorial image]]
Note that, with respect to the Hartree kernel (RPA calculations) the libxc is called to calculate the xc kernel within the ALDA, in particular the Perdew-Zunger parametrization for the correlation part of the LDA is used.  


Clearly there is very little difference between the two spectra. This highlights an important point in calculating excited state properties: generally, fewer G-vectors are needed than what are needed in DFT calculations. Regarding the spectrum itself, the first peak occurs at about 4.4eV. This is consistent with the minimum direct gap reported by Yambo: 4.28eV. The comparison with experiment (not shown) is very poor however.
Let's compare the absorption with and without the ALDA (+ local fields) included. By inspecting the o-3D-ALDA.eps_q1_inv_rpa_dyson file we find that this information is given in the 2nd and 4th columns, respectively:


If you make some mistake, and cannot reproduce this figure, you should check the value of <code>[[Variables#FFTGvecs|FFTGvecs]]</code> in the input file, delete the ''6Ry'' folder, and try again - taking care to plot the right file! (e.g. ''o-6Ry.eps_q1_ip_01'').
$ head -n30 o-3D_ALDA.eps_q1_inv_alda_dyson
# Absorption @ Q(1) [q->0 direction] : 0.7071068  0.7071068  0.0000000
#  E/ev[1]    EPS-Im[2]   EPS-Re[3]   EPSo-Im[4]  EPSo-Re[5]


==q-direction==
Plot the result:
Now let's select a different component of the dielectric tensor:
$ yambo -F yambo.in_IP -J 6Ry -V RL -o c
...
% [[Variables#LongDrXd|LongDrXd]]
'''0.000000''' | 0.000000 | '''1.000000''' |        # [Xd] [cc] Electric Field
%
...
$ yambo -F yambo.in_IP -J 6Ry -V RL
This time yambo reads from the ''6Ry'' folder, so it does not need to compute the dipole matrix elements again, and the calculation is fast. Plotting gives:
$ gnuplot
gnuplot> plot "o-6Ry.eps_q1_ip" t "q || x-axis" w l,"o-6Ry.eps_q1_ip_01" t "q || c-axis" w l


[[File:CH-hBN-ac.png|none|500px|Yambo tutorial image]]
$ gnuplot
The absorption is suppressed in the stacking direction. As the interplanar spacing is increased, we would eventually arrive at the absorption of the BN sheet (see [[Local fields]] tutorial).
gnuplot> plot "o-3D_ALDA.eps_q1_inv_alda_dyson" u 1:2 w l t 'ALDA',"o-3D_ALDA.eps_q1_inv_alda_dyson" u 1:4 w l t 'IPA'


==Non-local commutator==
[[File:3D_ALDA_hBN.png|none|500px|Yambo tutorial image]]
Last, we show the effect of switching off the non-local commutator term (see ''[Vnl,r]'' in the equation at the top of the page) due to the pseudopotential. As there is no option to do this inside yambo, you need to hide the database file. Change back to the ''q || (1 0 0)'' direction, and launch yambo with a different <code>-J</code> option:
$ mv SAVE/ns.kb_pp_pwscf SAVE/ns.kb_pp_pwscf_'''OFF'''
$ yambo -F yambo.in_IP -J '''6Ry_NoVnl''' -o c
$ yambo -F yambo.in_IP -J 6Ry_NoVnl
Note the warning in the output:
<---> [WARNING] Missing non-local pseudopotential contribution
which also appears in the report file, and noted in the database as <code>[r,Vnl] included      :no</code>. The difference is tiny:
[[File:CH-hBN-Vnl.png|none|500px|Yambo tutorial image]]


However, when your system is larger, with more projectors in the pseudopotential or more k-points (see the BSE tutorial), the inclusion of ''Vnl'' can make a huge difference in the computational load, so it's always worth checking to see if the terms are important in your system.
What is clear from this plot is that the ALDA kernel does not significantly change the optical spectrum when compared with the independent particle level of approximation (IPA). The same trend can be observed for other semiconductors and insulators.


==Summary==
==Summary==
From this tutorial you've learned:
From this tutorial you've learned:
* How to compute a simple optical spectrum  
* How to compute an optical spectrum within the adiabatic local density approximation
* How to reduce the computational load through reducing the G-vector/energy cut off and removing the Vnl term
* How to plot different components of the dielectric tensor
* How to use the <code>-J</code> option to neatly label and organise files and databases


== Links ==
== Links ==
* Next module: [[Local fields|Local fields]]
* Next module: [[Time-dependent density functional theory (long-range corrected kernel) | Long-range corrected kernel]]
* Back to [[First steps: a walk through from DFT to optical properties|First steps with yambo]] tutorial
* [[Tutorials|Back to tutorials menu]]
* [[Tutorials|Back to tutorials menu]]
== References ==

Latest revision as of 09:36, 7 November 2019

In this module you will learn how to calculate optical spectra at the time-dependent density-functional theory level for bulk hBN.

Background

Within the time-dependent DFT formalism the macroscopic dielectric function is obtained by solving the Dyson equation for the susceptibility X. In reciprocal space this is given by:

Yambo tutorial image

The f xc term is called the exchange-correlation (xc) kernel and describes the response of the xc potential at a time t to changes in the density at all previous times. Approximations for the xc kernel can be derived by taking the functional derivative of common approximations for the xc potential with respect to the density with the additional approximation that only instantaneous changes to the density are considered (adiabatic approximation).

Within this family we consider here the adiabatic local density approximation (ALDA)

Eq ALDA.png

Prerequisites

Cheatsheet on TDDFT
  • You must first complete the "How to use Yambo" and the "Optics at the independent particle level" and "local fields" modules

You will need:

  • The SAVE databases for bulk hBN
  • The yambo executable
  • gnuplot, for plotting spectra

Choosing input parameters

Enter the folder for bulk hBN that contains the SAVE directory, and generate the input file. From yambo -H you should understand that the correct option is yambo -o c -k <opt> . Where <opt> can be 'hartree' as in the tutorial on the local fields, or 'alda' or 'lrc' (see next tutorial) respectively for the ALDA and LRC approximations for the xc kernel outlined above. Let's start by generating the input for a calculation with the ALDA kernel: 

$ cd YAMBO_TUTORIALS/hBN-3D/YAMBO
$ yambo        (Initialization)
$ yambo -F yambo.in_ALDA -J 3D_ALDA -o c -k alda

We thus use a new input file yambo.in_ALDA, switch on the FFTGvecs variable, and label all outputs/databases with a 3D_ALDA tag. Make sure to set/modify all of the following variables to: 

Chimod= "ALDA"            # [X] IP/Hartree/ALDA/LRC/BSfxc
FxcGRLc=     3        Ry    # [Xd] Response block size
NGsBlkXd=     3        Ry    # [Xd] Response block size
% QpntsRXd
 1 |  1 |                   # [Xd] Transferred momenta
%
% EnRngeXd
 0.00000 | 10.00000 | eV    # [Xd] Energy range
%
% DmRngeXd
0.200000 | 0.200000 | eV    # [Xd] Damping range
%
ETStpsXd= 1001               # [Xd] Total Energy steps
% LongDrXd
1.000000 | 1.000000 | 0.000000 |        # [Xd] [cc] Electric Field
%

In this input file, we have selected:

  • the long-wavelength limit q = 0 (optical limit)
  • A q parallel to the BN planes
  • The ALDA kernel, and expanded G-vectors in the screening and in the xc kernel up to 3 Ry (about 85 G-vectors)

Optics runlevel

Save the input file and launch the code, keeping the command line options as before (i.e., just remove the lower case options): 

$ yambo -F yambo.in_ALDA -J 3D_ALDA
...
<01s> [05] Optics
<01s> [LA] SERIAL linear algebra
<01s> [WF] Performing Wave-Functions I/O from ./SAVE
<01s> [FFT-Rho] Mesh size:  9   9  21
<01s> [xc] Functional Slater exchange(X)+Perdew & Zunger(C)
<01s> [xc] LIBXC used to calculate xc functional
...

Note that, with respect to the Hartree kernel (RPA calculations) the libxc is called to calculate the xc kernel within the ALDA, in particular the Perdew-Zunger parametrization for the correlation part of the LDA is used.

Let's compare the absorption with and without the ALDA (+ local fields) included. By inspecting the o-3D-ALDA.eps_q1_inv_rpa_dyson file we find that this information is given in the 2nd and 4th columns, respectively:

$ head -n30 o-3D_ALDA.eps_q1_inv_alda_dyson

  1. Absorption @ Q(1) [q->0 direction] : 0.7071068 0.7071068 0.0000000
  2. E/ev[1] EPS-Im[2] EPS-Re[3] EPSo-Im[4] EPSo-Re[5]

Plot the result:

$ gnuplot gnuplot> plot "o-3D_ALDA.eps_q1_inv_alda_dyson" u 1:2 w l t 'ALDA',"o-3D_ALDA.eps_q1_inv_alda_dyson" u 1:4 w l t 'IPA'

Yambo tutorial image

What is clear from this plot is that the ALDA kernel does not significantly change the optical spectrum when compared with the independent particle level of approximation (IPA). The same trend can be observed for other semiconductors and insulators.

Summary

From this tutorial you've learned:

  • How to compute an optical spectrum within the adiabatic local density approximation

Links

Retrieved from "https://wiki.yambo-code.eu/wiki/index.php?title=Time-dependent_density_functional_theory_(standard_kernel)&oldid=2910"