Optics at the independent particle level: Difference between revisions

From The Yambo Project
Jump to navigation Jump to search
 
(35 intermediate revisions by 3 users not shown)
Line 1: Line 1:
In this tutorial you will learn how to calculate optical spectra at the RPA or independent particle level for bulk hBN.
In this tutorial you will learn how to calculate optical spectra at the independent particle level for bulk hBN.
 
== Background ==
[[File:CH1.png|x50px|Yambo tutorial image]]
The dielectric function in the long-wavelength limit, at the independent particle level (RPA without local fields), is essentially given by the following:
[[File:CH1.png|none|x50px|Yambo tutorial image]]
In practice, Yambo does not use this expression directly but solves the Dyson equation for the susceptibility ''X'', which is described in the [[Local fields]] module.


== Prerequisites ==
== Prerequisites ==
 
[[File:Yambo-handbook-v4.1.2-p-3.png|thumb|Cheatsheet on linear response|150px]]
* You must first complete the "How to use Yambo" tutorial
* You must first complete the [[Modules#How to use Yambo|How to use Yambo]] modules


'''You will need''':
'''You will need''':
Line 13: Line 15:


== Choosing input parameters ==
== Choosing input parameters ==
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</code>. Let's add some command line options:
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 -optics c</code> (or <code>yambo -o c</code>). Let's add some command line options:
<!--If you have already worked in this folder, you might like to clean up any old files, and run the initialization again:
<!--If you have already worked in this folder, you might like to clean up any old files, and run the initialization again:
  cd YAMBO_TUTORIALS/hBN/YAMBO
  cd YAMBO_TUTORIALS/hBN/YAMBO
Line 20: Line 23:
Take care to '''not''' delete the <code>ns.*</code> files. -->
Take care to '''not''' delete the <code>ns.*</code> files. -->
  $ cd YAMBO_TUTORIALS/hBN/YAMBO
  $ cd YAMBO_TUTORIALS/hBN/YAMBO
  $ yambo -o c -F yambo.in_IP -J Full
  $ yambo               ''(initialization)''
This corresponds to optical properties in G-space at the independent particle level (<code>Chimod= "IP"</code>).
$ yambo -F yambo.in_IP -o c
====Optics runlevel====
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>).
Let's calculate just for the long-wavelength limit ''q = 0'' which is the valid case for optical spectra. This always corresponds to the ''first'' q-point. Change the following variables in the input file to:
 
==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]]
  % [[Variables#QpntsRX|QpntsRXd]]
   1 |  '''1''' |                  # [Xd] Transferred momenta
   1 |  '''1''' |                  # [Xd] Transferred momenta
Line 32: Line 38:
  $ yambo -F yambo.in_IP -J Full
  $ yambo -F yambo.in_IP -J Full
  ...
  ...
  <---> [05] Optics
  <---> [04] Dipoles
  <---> [LA] SERIAL linear algebra
  <---> [DIP] Checking dipoles header
  <---> [DIP] Checking dipoles header
<---> [WARNING] DIPOLES database not correct or missing
  <---> [x,Vnl] computed using 4 projectors
  <---> [x,Vnl] computed using 4 projectors
  <---> [M  0.017 Gb] Alloc WF ( 0.016)
  <04s> Dipoles: P, V and iR (T) |########################################| [100%] 04s(E) 04s(X)
<---> [WF] Performing Wave-Functions I/O from ./SAVE
  <04s> [DIP] Writing dipoles header
<01s> Dipoles: P and iR (T): |########################################| [100%] 01s(E) 01s(X)
  <04s> [WARNING] Only commutator with non-local pseudopotential contribution included in position and velocity dipol
  <01s> [M 0.001 Gb] Free WF ( 0.016)
<04s> [05] Optics
  <01s> [DIP] Writing dipoles header
<04s> [LA] SERIAL linear algebra
  <01s> [X-CG] R(p) Tot o/o(of R) :   5501   52992    100
  <04s> [DIP] Checking dipoles header
  <01s> Xo@q[1] |########################################| [100%] --(E) --(X)
  <04s> [X-CG] R(p) Tot o/o(of R):   5499   52992    100
  <01s> [06] Game Over & Game summary
  <05s> Xo@q[1] |########################################| [100%] --(E) --(X)
  <05s> [06] Timing Overview
<05s> [07] Memory Overview
<05s> [08] Game Over & Game summary
 
  $ ls
  $ ls
  Full  SAVE  yambo.in_IP  r_setup
  Full  SAVE  yambo.in_IP  r_setup
  o-Full.eel_q1_ip  o-Full.eps_q1_ip  r-Full_optics_chi
  o-Full.eel_q1_ip  o-Full.eps_q1_ip  r-Full_optics_*
Let's take a moment to understand what Yambo has done in side the Optics runlevel:
Let's take a moment to understand what Yambo has done inside the Optics runlevel:
# Compute the ''[x,Vnl]'' term
# Compute the ''[x,Vnl]'' term
# Read the wavefunctions from disc [WF]
# Read the wavefunctions from disc [WF]
# Compute the ''dipoles'', i.e. matrix elements of '''p'''
# Compute the ''dipoles'', i.e. matrix elements of '''p'''
# Write the dipoles to disk as <code>SAVE/ndb.dip*</code> databases. This you can see in the report file:
# 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
  $ grep -A20 "WR" r-Full_optics_dipoles_chi
yambo 4.5 (or before)
  [WR./Full//ndb.dip_iR_and_P]
  [WR./Full//ndb.dip_iR_and_P]
  Brillouin Zone Q/K grids (IBZ/BZ):  14  72  14  72
  Brillouin Zone Q/K grids (IBZ/BZ):  14  72  14  72
Line 62: Line 73:
  [r,Vnl] included      :yes
  [r,Vnl] included      :yes
  ...
  ...
<ol start="4">
yambo 5.0
<li>Finally, Yambo computes ''X<sub>0</sub>'' for this ''q'', and writes the dielectric function inside the <code>o-Full.eps_q1_ip</code> file for plotting</li>
[WR./Full//ndb.dipoles]---------------------------------------------------------
  Brillouin Zone Q/K grids (IBZ/BZ)                :  14  72  14  72
  RL vectors                                      :  1491 [WF]
  Fragmentation                                    : yes
  Electronic Temperature                          :  0.000000 [K]
  Bosonic    Temperature                          :  0.000000 [K]
  DIP band range                                  :    1  100
  DIP band range limits                            :  8  9
  DIP e/h energy range                            : -1.000000 -1.000000 [eV]
  RL vectors in the sum                            :  1491
  [r,Vnl] included                                : yes
  Bands ordered                                    : yes
  Direct v evaluation                              : no
  Approach used                                    : G-space v
  Dipoles computed                                : R V P
  Wavefunctions                                    : Slater exchange(X)+Perdew & Zunger(C)
- S/N 003471 ---------------------------------------------- v.05.01.01 r.21242 -
<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>
</ol>


====Cut off====
==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:
Before plotting the output, let's change a few more variables. The previous calculation used ''all'' the G-vectors in expanding the wavefunctions, up to 1491 (~1016 components). 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
  $ yambo -F yambo.in_IP '''-J 6Ry -V RL''' -o c
Change the '''value''' of <code>[[Variables#FFTGvecs|FFTGvecs]]</code> and also its '''unit''':
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
  [[Variables#FFTGvecs|FFTGvecs]]= '''6'''          '''Ry'''    # [FFT] Plane-waves
Save the input file and launch the code again:
Save the input file and launch the code again:
   $ yambo -F yambo.in_IP '''-J 6Ry -V RL'''
   $ yambo -F yambo.in_IP '''-J 6Ry'''
and then plot the ''o-Full.eps_q1_ip'' and ''o-6Ry.eps_q1_ip'' files:
and then plot the ''o-Full.eps_q1_ip'' and ''o-6Ry.eps_q1_ip'' files:
  $ gnuplot
  $ gnuplot
  gnuplot> plot "o-Full.eps_q1_ip" w l,"o-6Ry.eps_q1_ip" w p
  gnuplot> plot "o-Full.eps_q1_ip" w l,"o-6Ry.eps_q1_ip" w p


[[File:CH-hBN-6Ry.png|500px|Yambo tutorial image]]
[[File:CH-hBN-6Ry.png|none|500px|Yambo tutorial image]]
 
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.


Clearly there is very little difference between the two spectra.
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'').


====q-direction====
==q-direction==
Now let's select a different component of the dielectric tensor:
Now let's select a different component of the dielectric tensor:
  $ yambo -o c -F yambo.in_IP -J 6Ry -V RL
  $ yambo -F yambo.in_IP -J 6Ry -V RL -o c
  ...
  ...
  % LongDrXd
  % [[Variables#LongDrXd|LongDrXd]]
  '''0.000000''' | 0.000000 | '''1.000000''' |        # [Xd] [cc] Electric Field
  '''0.000000''' | 0.000000 | '''1.000000''' |        # [Xd] [cc] Electric Field
  %
  %
  ...
  ...
  $ yambo -F yambo.in_IP -J 6Ry -V RL
  $ yambo -F yambo.in_IP -J 6Ry
This time yambo reads from the 6Ry folder, so it does not need to compute the dipoles again, and the calculation is fast. Plotting gives:
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
  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
  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|500px|Yambo tutorial image]]
[[File:CH-hBN-ac.png|none|500px|Yambo tutorial image]]
====Non-local commutator====
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).
Last, we show the effect of switching off the non-local commutator term (see the expression at the top) due to the pseudopotential. As there is no option to do this inside yambo, you need to rename the database file. Change back to the ''q || (1 0 0)'' direction, and launch yambo with a different <code>-J</code> option:
 
==Non-local commutator==
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'''
  $ mv SAVE/ns.kb_pp_pwscf SAVE/ns.kb_pp_pwscf_'''OFF'''
  $ yambo -o c -F yambo.in_IP -J '''6Ry_NoVnl'''
  $ yambo -F yambo.in_IP -J '''6Ry_NoVnl''' -o c
change back
% LongDrXd
1.000000 | 0.000000 | 0.000000 |        # [Xd] [cc] Electric Field
%
and run
  $ yambo -F yambo.in_IP -J 6Ry_NoVnl
  $ yambo -F yambo.in_IP -J 6Ry_NoVnl
Note the warning in the output:
Note the warning in the output:
  <---> [WARNING] Missing non-local pseudopotential contribution
  <---> [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:
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|500px|Yambo tutorial image]]
[[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.
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.


====Summary====
==Summary==
From this tutorial you've learned:
From this tutorial you've learned:
* How to compute a simple optical spectrum  
* How to compute a simple optical spectrum  
Line 115: Line 153:
* How to use the <code>-J</code> option to neatly label and organise files and databases  
* How to use the <code>-J</code> option to neatly label and organise files and databases  


== Links ==
<br>
* Next module: [[Local fields|Local fields]]
{| style="width:100%" border="1"
* [[Tutorials|Back to tutorials menu]]
|style="width:15%; text-align:left"|Prev: Back to [[First_steps:_a_walk_through_from_DFT_to_optical_properties#Step_4:_Optical_absorption_in_hBN|First steps]]
|style="width:70%; text-align:center"|Now: [[Tutorials|Tutorials Home]] --> [[First_steps:_a_walk_through_from_DFT_to_optical_properties|First steps]] --> [[Bulk_material:_h-BN|bulk hBN]]  --> [[Optics_at_the_independent_particle_level|IP optics]]
|style="width:15%; text-align:right"|Next: Back to [[First_steps:_a_walk_through_from_DFT_to_optical_properties#Step_4:_Optical_absorption_in_hBN|First steps]]
|-
|}

Latest revision as of 09:26, 9 February 2023

In this tutorial you will learn how to calculate optical spectra at the independent particle level for bulk hBN.

Background

The dielectric function in the long-wavelength limit, at the independent particle level (RPA without local fields), is essentially given by the following:

Yambo tutorial image

In practice, Yambo does not use this expression directly but solves the Dyson equation for the susceptibility X, which is described in the Local fields module.

Prerequisites

Cheatsheet on linear response

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, run the initialization and generate the input file. From yambo -h you should understand that the correct option is yambo -optics c (or yambo -o c). Let's add some command line options:

$ cd YAMBO_TUTORIALS/hBN/YAMBO
$ yambo               (initialization)
$ yambo -F yambo.in_IP -o c

This corresponds to optical properties in G-space at the independent particle level: in the input file this is indicated by (Chimod= "IP").

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:

% QpntsRXd
 1 |  1 |                   # [Xd] Transferred momenta
%
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):

$ yambo -F yambo.in_IP -J Full
...
<---> [04] Dipoles
<---> [DIP] Checking dipoles header
<---> [WARNING] DIPOLES database not correct or missing
<---> [x,Vnl] computed using 4 projectors
<04s> Dipoles: P, V and iR (T) |########################################| [100%] 04s(E) 04s(X)
<04s> [DIP] Writing dipoles header
<04s> [WARNING] Only commutator with non-local pseudopotential contribution included in position and velocity dipol
<04s> [05] Optics
<04s> [LA] SERIAL linear algebra
<04s> [DIP] Checking dipoles header
<04s> [X-CG] R(p) Tot o/o(of R):    5499   52992     100
<05s> Xo@q[1] |########################################| [100%] --(E) --(X)
<05s> [06] Timing Overview
<05s> [07] Memory Overview
<05s> [08] 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_*

Let's take a moment to understand what Yambo has done inside the Optics runlevel:

  1. Compute the [x,Vnl] term
  2. Read the wavefunctions from disc [WF]
  3. Compute the dipoles, i.e. matrix elements of p
  4. Write the dipoles to disk as Full/ndb.dip* databases. This you can see in the report file:
$ grep -A20 "WR" r-Full_optics_dipoles_chi 

yambo 4.5 (or before)

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

yambo 5.0

[WR./Full//ndb.dipoles]---------------------------------------------------------
 Brillouin Zone Q/K grids (IBZ/BZ)                :   14   72   14   72
 RL vectors                                       :  1491 [WF]
 Fragmentation                                    : yes
 Electronic Temperature                           :  0.000000 [K]
 Bosonic    Temperature                           :  0.000000 [K]
 DIP band range                                   :    1  100
 DIP band range limits                            :   8   9
 DIP e/h energy range                             : -1.000000 -1.000000 [eV]
 RL vectors in the sum                            :  1491
 [r,Vnl] included                                 : yes
 Bands ordered                                    : yes
 Direct v evaluation                              : no
 Approach used                                    : G-space v
 Dipoles computed                                 : R V P
 Wavefunctions                                    : Slater exchange(X)+Perdew & Zunger(C)
- S/N 003471 ---------------------------------------------- v.05.01.01 r.21242 -
  1. Finally, Yambo computes the non-interacting susceptibility X0 for this q, and writes the dielectric function inside the o-Full.eps_q1_ip file for plotting

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, up to 1491 (~1016 components). 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 FFTGvecs and also its unit from RL (number of G-vectors) to Ry (energy in Rydberg):

FFTGvecs= 6           Ry    # [FFT] Plane-waves

Save the input file and launch the code again:

 $ yambo -F yambo.in_IP -J 6Ry

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
Yambo tutorial image

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.

If you make some mistake, and cannot reproduce this figure, you should check the value of FFTGvecs 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).

q-direction

Now let's select a different component of the dielectric tensor:

$ yambo -F yambo.in_IP -J 6Ry -V RL -o c
...
% LongDrXd
0.000000 | 0.000000 | 1.000000 |        # [Xd] [cc] Electric Field
%
...
$ yambo -F yambo.in_IP -J 6Ry

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
Yambo tutorial image

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

Non-local commutator

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 -J option:

$ mv SAVE/ns.kb_pp_pwscf SAVE/ns.kb_pp_pwscf_OFF
$ yambo -F yambo.in_IP -J 6Ry_NoVnl -o c

change back

% LongDrXd
1.000000 | 0.000000 | 0.000000 |        # [Xd] [cc] Electric Field
%

and run

$ 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 [r,Vnl] included :no. The difference is tiny:

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.

Summary

From this tutorial you've learned:

  • How to compute a simple optical spectrum
  • 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 -J option to neatly label and organise files and databases


Prev: Back to First steps Now: Tutorials Home --> First steps --> bulk hBN --> IP optics Next: Back to First steps