XANES in Diamond: Difference between revisions

From VASP Wiki
VisualWikitext
 
(33 intermediate revisions by 2 users not shown)
Line 1: Line 1:
(UNDER CONSTRUCTION)
Important: This feature will be only available from VASP 6.0.


{{Template:XAS}}
{{Template:XAS - Tutorial}}


== Task ==
== Task ==
Line 19: Line 19:
  0.25 0.25 0.25
  0.25 0.25 0.25


*We will not use this structure as input but rather use it to construct a super cell from it that is actually used in the calculations.
Above, we show the {{FILE|POSCAR}} file for the primitive unit cell of diamond. Note that we will not use this
structure as input for the calculation. Instead, we use it to construct a {{FILE|POSCAR}} file for a
<math>3 \times 3 \times 3</math> super cell actually used in the calculation.


=== {{FILE|INCAR}} ===
=== {{FILE|INCAR}} ===
Line 35: Line 37:
  {{TAGBL|LREAL}} = A
  {{TAGBL|LREAL}} = A


*To promote a core electron into the conduction bands and hence create the core-hole {{TAG|ICORELEVEL}}=2 has to be set. This corresponds to the final state approximation
*To promote a core electron into the conduction bands and hence create the core-hole {{TAG|ICORELEVEL}}=2 has to be set. This corresponds to the final state approximation.
*{{TAG|CLNT}}=1 selects the first atom species in the {{TAG|POSCAR}} file.
*{{TAG|CLNT}}=1 selects the first atom species in the {{TAG|POSCAR}} file.
*{{TAG|CLN}}=1 selects main quantum nuhttps://cms.mpi.univie.ac.at/wiki/index.mber 1 (hence K-edge).
*{{TAG|CLN}}=1 selects main quantum number 1 (hence K-edge).
*{{TAG|CLL}}=0 selects angular quantum number 0 (s).
*{{TAG|CLL}}=0 selects angular quantum number 0 (s).
*{{TAG|CLZ}}=1.0 selects the charge of the core hole. By setting this number to a fractional value we can mimick different screening of the electrons. Since this purely exploits error cancellation and the physical background of non-integer charges is not defined well, it should be only used with caution.
*{{TAG|CLZ}}=1.0 selects the charge of the core hole. By setting this number to a fractional value we can mimic different screenings of the electrons. These non-integer values should only be used with caution, since this purely exploits error cancellation and does not correspond to a physically correct description of the screening.
*By setting {{TAG|CH_LSPEC}}=''.TRUE.'' we enable the calculation of matrix elements between core and conduction states and the calculation of the core electron absorption spectrum.
*By setting {{TAG|CH_LSPEC}}=''.TRUE.'', we enable the calculation of matrix elements between core and conduction states and the calculation of the core electron absorption spectrum.
*The broadening of the core electron absorption spectrum is controlled by the tag {{TAG|CH_SIGMA}}. Usually it is good practice to set this value low and broaden the spectrum in post processing.
*The broadening of the core electron absorption spectrum is controlled by the tag {{TAG|CH_SIGMA}}. Usually it is good practice to set this value low and broaden the spectrum in post processing.
*We have to set {{TAG|NBANDS}} to a larger value to consider enough conduction band states in the calculation.
*We have to set {{TAG|NBANDS}} to a larger value to consider enough conduction band states in the calculation.
*Since super cells are used the calculation of the projection operators in real space is much faster, hence {{TAG|LREAL}}=''A'' is set.
*Since super cells are used the calculation of the projection operators in real space ({{TAG|LREAL}}=''A'') is much faster.


== Calculation ==
== Calculation ==
Line 49: Line 51:
=== Step 1 build a supercell ===
=== Step 1 build a supercell ===


To perform calculations with a core~/VASP_WIKI/WIKI_TUTORIALS/CORE_HOLES/Diamond/S3x3x3/FCH/P4V-hole we need sufficiently large super cells to reduce the interaction of neighbouring core-holes. Usually the obtained spectrum with respect to the size of the super cell has to be converged. This means calculations are done with successively larger cells until no significant change is visible in the spectrum anymore. To save computational time we chose a <math>3\times3\times3</math> cell for this tutorial, although for converged values one should use at least <math>4\times4\times4</math>.  
In the periodic boundary conditions of VASP, the core-hole interacts with its periodic replica so that we need sufficiently large super cells to reduce this spurious interaction. To this end, we employ successively large cells until the spectrum shows no significant changes. For this tutorial, we illustrate the calculation of a core-hole using a <math>3\times3\times3</math> cell to allow for a reasonably fast calculation. However, for converged values one should use at least <math>4\times4\times4</math> super cell.  


The super cell can be obtained either by taking the file POSCAR.3x3x3 provided with this tutorial. In the following we show how to get the super cell from the primitive cell using p4vasp:
The super cell can be obtained by either taking the file POSCAR.3x3x3 provided with this tutorial or constructing the {{FILE|POSCAR}} file from the primitive cell using p4vasp, which is demonstrated below:
*Open p4vasp by typing ''p4v'' on the terminal.
*Open p4vasp by typing ''p4v'' on the terminal.
*Load the primitive cell by clicking on '''File'''&rarr;'''Load system''':
*Load the primitive cell by clicking on '''File'''&rarr;'''Load system''':
Line 61: Line 63:


=== Step 2 Prepare input files ===
=== Step 2 Prepare input files ===
The first few lines of the {{TAG|POSCAR}} file for the super cell should look like the following
The first few lines of the generated {{TAG|POSCAR}} file for the super cell should look like the following
  cubic diamond
  cubic diamond
  3.567
  3.567
   +1.5000000000  +1.5000000000  +0.0000000000
   +1.5000000000  +1.5000000000  +0.0000000000
   +0.0000000000  +1.5000000000  +1.~/VASP_WIKI/WIKI_TUTORIALS/CORE_HOLES/Diamond/S3x3x3/FCH/P4V5000000000
   +0.0000000000  +1.5000000000  +1.5000000000
   +1.5000000000  +0.0000000000  +1.5000000000
   +1.5000000000  +0.0000000000  +1.5000000000
   54
   54
Line 73: Line 75:
   +0.5000000000  +0.0000000000  +0.5000000000
   +0.5000000000  +0.0000000000  +0.5000000000
   ...
   ...
Here we have only one atom species with 54 atoms (line 6). To construct a core hole on a single atom we have to, take one of these atoms and treat it as a different species. We choose the first atom and we end up with 1 and 53 in the 6th line of the modified {{TAG|POSCAR}} file
Here, all the 54 atoms are of the same species (line 6). To distinguish between the atom with the core-hole and the rest, we treat one atom as a different species. Choosing the first atom, we replace the 54 in the 6th line with 1 and 53 and obtain the following {{TAG|POSCAR}} file
  cubic diamond
  cubic diamond
  3.567
  3.567
Line 85: Line 87:
   +0.5000000000  +0.0000000000  +0.5000000000
   +0.5000000000  +0.0000000000  +0.5000000000
   ...
   ...
Accordingly we have to set {{TAGBL|CLNT}}=1 in the {{TAG|INCAR}} file which selects the first atom species in {{TAG|POSCAR}} to carry the core-hole. Also one has to create a {{TAG|POTCAR}} file with the PAW/PS information for both species. Since the two species are both diamond this can be simply done by concatenation of the {{TAG|POTCAR}} file for diamond by using the command ''cat POT_C POT_C > POTCAR''. Alternatively the {{TAG|POTCAR}} appropriate for core-hole calculations in diamond is provided in the tar file.
In the {{TAG|INCAR}} file, we specify that the first species carries the core-hole by setting {{TAGBL|CLNT}}=1. We create a {{TAG|POTCAR}} file with the PAW/PS information for both species. Since both species are carbon this amounts simply to the concatenation of the {{TAG|POTCAR}} file for carbon
  cat POT_C POT_C > POTCAR
We provide the resulting {{TAG|POTCAR}} file for this core-hole calculation in the tar file of this tutorial.


The number of used bands in the calculation have to be set manually. Usually one needs to select enough bands, depending on how far one wants to calculate the spectrum and on the number of electrons in the system, to have enough conduction states available in the calculation. Selecting too high numbers a priori is also not good, since it drastically increases the computation time. Hence want has to test the optimal number of bands. In our example we set {{TAG|NBANDS}}=300.
To calculate accurate spectra, we need to include a sufficient number of conduction states. The required number of bands depends on the number of electrons in the system and the energy range of the spectrum. We can manually adjust the number of bands with the {{TAG|NBANDS}} variable. However, since the computation time increases drastically with the number of bands, selecting initially very large numbers is also not advisable. Hence, one has to increase the number of bands to find the optimum of computational effort and accuracy. In this tutorial, we use {{TAG|NBANDS}}=300.


The rest of the {{TAG|INCAR}} file is described above.
The other input variables in the {{TAG|INCAR}} file are described above.


'''Mind''': The multiplicity of the species carrying the core hole has to be 1 otherwise the code will not work correctly. Also mind that the selected species ({{TAGBL|CLNT}} in the {{TAG|INCAR}} file) is consistent with the order of the species specified in the {{TAG|POSCAR}} and {{TAG|POTCAR}} files.
'''Mind''': The multiplicity of the species carrying the core-hole has to be 1 otherwise the code will not work. Also mind that the selected species ({{TAGBL|CLNT}} in the {{TAG|INCAR}} file) is consistent with the order of the species specified in the {{TAG|POSCAR}} and {{TAG|POTCAR}} files.




=== Step 3 Running Calculation ===
=== Step 3 Running Calculation ===


The SCF calculation with the core-hole and afterwards the calculation dielectric matrix (spectrum) is done in a single calculation.
Both the SCF calculation with the core-hole and the subsequent calculation of the dielectric matrix (spectrum) are done in a run of VASP. To minimize the spurious interaction between core-holes in neighboring cells requires large super cells and to reduce the computational time it is advisable to run a parallel VASP calculation
Usually the command to run the calculation looks like this: ''mpirun -np $np vasp_version'', where ''$np'' corresponds to the number of processes and ''_version'' in the executable usually stands for std, gam, ncl namely standard, gamma-point only and non-collinear version, respectively.  
  mpirun -np $np vasp_version
Since we have to use super cells in the calculation, which are large enough to minimize the interaction between core-holes on neighbouring cells, parallel execution with many computational cores is necessary to achieve good results within reasonable time. Hence for this example we use a <math>3\times3\times3</math> cell, which gives sufficiently accurately results but more importantly finishes fast using only a few number of cores. For accurate results one should use at least a <math>4\times4\times4</math> cell.
Here, ''$np'' corresponds to the number of processes and the ''_version'' in the executable stands for ''std'', ''gam'', and ''ncl'': the standard, <math>\Gamma</math>-point only, and non-collinear version, respectively.
The <math>3\times3\times3</math> cell used in this tutorial gives qualitatively correct results and can be completed even with a small number of processes. You can verify the spectrum on the larger <math>4\times4\times4</math> cell, which we provide in the tar file of this tutorial.


=== Step 4 Extraction of XAS Spectrum ===
=== Step 4 Extraction of XAS Spectrum ===


To get the spectrum we want to plot the imaginary part of the frequency depedent dielectric function. This is written out in the {{TAG|OUTCAR}} file
The XAS spectrum is proportional to the imaginary part of the frequency-dependent dielectric function, which is written in the {{TAG|OUTCAR}} file
   frequency dependent IMAGINARY DIELECTRIC FUNCTION (independent particle, no local field effects) density-density
   frequency dependent IMAGINARY DIELECTRIC FUNCTION (independent particle, no local field effects) density-density
       E(ev)      X        Y        Z        XY        YZ        ZX
       E(ev)      X        Y        Z        XY        YZ        ZX
Line 111: Line 116:
   ...
   ...


Usually we are interested in the sum of all components of the dielectric matrix. You can either obtain this by your script or you use the script provided in this example. If you decide for the latter use the command ''perl ./plot_core_imdiel.pl''. This will create the file CORE_DIELECTRIC_IMAG.dat which contains the summed up imaginary part of the dielectric matrix.
Usually we are interested in the sum of all components of the dielectric matrix. You can obtain this by the script provided in this tutorial ''plot_core_imdiel.sh''
 
<div class="toccolours mw-customtoggle-script">'''Click to show/''plot_core_imdiel.sh''</div>
<div class="mw-collapsible mw-collapsed" id="mw-customcollapsible-script">
<pre>
#!/bin/bash
 
parallel=-1
normal=-1
all=-1
tauc=-1
trace=-1
while [[ $# -gt 0 ]]
do
  key="$1"
  case $key in
      -parallel) parallel=0
      ;;
      -normal) normal=0
      ;;
      -trace) trace=0
      ;;
      -tauc) tauc=0
      ;;
  esac
  shift
done
 
 
cat > helpscript.perl  <<EOF
#!/bin/perl
 
use strict;
use warnings;
my \$mode=shift;
 
while(<>)
{
  chomp;
  if(\$_ =~ /frequency dependent IMAGINARY DIELECTRIC FUNCTION/)
  {
      \$_=<>;
      \$_=<>;
      while (<>)
      {
        my \$sum=0;
        if (\$_ !~ /[0-9]/) {last;}
        chomp;
        \$_=~s/^/ /;
        my @help=split(/[\t,\s]+/);
        if (\$help[2]=~/NaN/||\$help[3]=~/NaN/||\$help[4]=~/NaN/) {next;}
        if (\$help[5]=~/NaN/||\$help[6]=~/NaN/||\$help[4]=~/NaN/) {next;}
        if (\$mode==0) {\$sum=\$help[2]+\$help[3]+\$help[4]+\$help[5]+\$help[6]+\$help[7];}
        if (\$mode==1) {\$sum=\$help[4];}
        if (\$mode==2) {\$sum=\$help[2]+\$help[3];}
        if (\$mode==3) {\$sum=\$help[2]+\$help[3]+\$help[4];}
        if (\$mode==4) {\$sum=(\$help[1]*\$help[1]*(\$help[2]+\$help[3]+\$help[4]+\$help[5]+\$help[6]+\$help[7]))**0.5;}
        if (\$mode==5) {\$sum=(\$help[1]*\$help[1]*(\$help[4]))**0.5;}
        if (\$mode==6) {\$sum=(\$help[1]*\$help[1]*(\$help[2]+\$help[3]))**0.5;}
        if (\$mode==7) {\$sum=(\$help[1]*\$help[1]*(\$help[2]+\$help[3]+\$help[4]))**0.5;}
        print \$help[1]," ",\$sum,"\n";
      }
  }
  last if eof;
}
EOF
 
if [[ $normal -eq 0 ]]; then
  if [[ $tauc -eq 0 ]]; then
      perl helpscript.perl 4 OUTCAR > CORE_DIELECTRIC_IMAG.dat
  else
      perl helpscript.perl 1 OUTCAR > CORE_DIELECTRIC_IMAG.dat
  fi
else
  if [[ $parallel -eq 0 ]]; then
      if [[ $tauc -eq 0 ]]; then
        perl helpscript.perl 5 OUTCAR > CORE_DIELECTRIC_IMAG.dat
      else
        perl helpscript.perl 2 OUTCAR > CORE_DIELECTRIC_IMAG.dat
      fi
  else
      if [[ $trace -eq 0 ]]; then
        if [[ $tauc -eq 0 ]]; then
            perl helpscript.perl 6 OUTCAR > CORE_DIELECTRIC_IMAG.dat
        else
            perl helpscript.perl 3 OUTCAR > CORE_DIELECTRIC_IMAG.dat
        fi
      else
        if [[ $tauc -eq 0 ]]; then
            perl helpscript.perl 7 OUTCAR > CORE_DIELECTRIC_IMAG.dat
        else
            perl helpscript.perl 0 OUTCAR > CORE_DIELECTRIC_IMAG.dat
        fi
      fi
  fi
fi
rm helpscript.perl
 
</pre>
</div>


For this example we want to compare with experimental and theoretical XAS spectra from literature. They are provided in the files ''C_XAS_aligned_to_VASP.dat'' and ''C_PARATEC_aligned_to_VASP.dat''. The literature calculations were obtained using the PARATEC code, which is also a PAW/Pseudopotential code. Both spectra are shifted to coincide with the VASP at the first peak. One can either plot these two spectra together with the VASP spectrum using his own prefered method of choice or using the script provided with this tutorial and invoking the command ''gnuplot gnuplot.script''. The output should hopefully look like the this:
To use it type:
[[File:Fig XAS 4.png]]
bash ./plot_core_imdiel.sh


The experimental offset of the spectrum (absolute value of the first peak) is rather impossible to reproduce, since the calculated core energies are very different compared to experiment. Usually there is even a noticable deviation between calculations using different codes. It is accepted in literature to look at the relative peak positions in the spectra and the spectra can be scaled arbitrarily . In this example we scaled the experiment to VASP, since in this way the obtained results can be very easily compared to the experiment using a script. Usually one would either scale the first peak to 0 or would scale the calculated value to the experiment. Additionally the intensity of the spectrum can be scaled arbitrarily. So in this example we align the position and the height of the first peak for the calculations. The experiment is a
This will create the file CORE_DIELECTRIC_IMAG.dat containing the sum of the imaginary part of the dielectric matrix.
Note that the absolute values of the experimental peak positions is not captured due to fundamental limitations of local DFT to describe the core level energies. Usually, there is even a noticeable deviation between calculations using different codes. Therefore, it is accepted in literature to look at the relative peak positions in the spectra and their relative intensity.
 
We compare the results obtained with VASP to experimental<ref name="Ma.PRL"/> and theoretical<ref name="Tallefumier.PRB"/> XAS spectra from literature provided in the files ''C_XAS_aligned_to_VASP.dat'' and ''C_PARATEC_aligned_to_VASP.dat'', respectively. The theoretical reference calculation was obtained using the PARATEC code and relies on PAW/Pseudopotential similar to VASP. We provide a script with this tutorial to compare these literature results to the spectrum obtained with VASP:
 
<div class="toccolours mw-customtoggle-script">'''Click to show/''gnuplot_XANES_C.script''</div>
<div class="mw-collapsible mw-collapsed" id="mw-customcollapsible-script">
<pre>
unset ytics
set xrange [280:310]
set xlabel "Energy (eV)"
set ylabel "Absorption (arbitrary units)"
plot "CORE_DIELECTRIC_IMAG.dat" using 1:2 with lines lw 2 ti "VASP",\
    "C_PARATEC_aligned_to_VASP.dat" using 1:2 with lines ti "PAW lit",\
    "C_XAS_aligned_to_VASP.dat" using 1:2 with lines ti "Exp"
pause -1
</pre>
</div>
 
To use that script type:
  gnuplot gnuplot_XANES_C.script
The file ''plot.sh'' constitutes a convenient wrapper around these post processing steps. The resulting spectra should look like this:<br />
[[File:Fig XAS 4.png|600px]]
 
Because DFT cannot reproduce the absolute position (see above), we have shifted both spectra so that the position of the first peak  coincides with the VASP. In this example, we scaled the experiment to VASP, since in this way the obtained results can be very easily compared using a script. Usually one would either scale the first peak to 0 or would scale the calculated value to the experiment. Additionally the intensity of the spectrum can be scaled arbitrarily. So in this example, we align the position and the height of the first peak for the calculations. The experiment is a little bit more tricky. It's a matter of taste what to consider as the first peak, but we decided that most likely the second peak corresponds to the first peak in the calculations and the first peak in experiment is a shoulder that is simply not pronounced in the calculations.
 
Another important issue is the broadening. Because the observed broadening is driven by many factors depending on the particular experimental setup, it is not possible to reproduce the broadening exactly. Therefore, we arbitrarily choose a broadening that gives approximately the same width as experiment. For simplicity in this tutorial, we choose to use a 0.5 eV constant Gaussian broadening (the broadening is determined by the {{TAG|ISMEAR}} tag and we used {{TAG|ISMEAR}}=0 which corresponds to a Gaussian broadening). For more elaborate spectra we strongly advise users to choose a 0.05 eV broadening and apply the desired broadening as post-processing.
 
Apart from the lower broadening width, we get a quite reasonable agreement with the theoretical literature calculation. We stress again that the <math>3\times3\times3</math> super cell in this example are not fully converged. The interested user can repeat the calculations for a larger <math>4\times4\times4</math> super cell. The files for this example are also given in the tar file. Be aware that the larger number of atoms in the bigger super cell requires an adjustment of the {{TAG|NBANDS}} variable..


== Download ==
== Download ==
[[Media:XANES_Diamond.tgz| XANES_Diamond.tgz]]
== References ==
<references>
<ref name="Ma.PRL">[https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.69.2598 Y.Ma et al., Phys. Rev. Lett 69, 2598 (1992).]</ref>
<ref name="Tallefumier.PRB">[https://journals.aps.org/prb/abstract/10.1103/PhysRevB.66.195107 M.Tallefumier et al., Phys. Rev. B 66, 195107 (2002).]</ref>
</references>


{{Template:XAS}}
{{Template:XAS}}
Back to the [[The_VASP_Manual|main page]].


[[Category:Examples]]
[[Category:Examples]]

Latest revision as of 15:30, 15 April 2022

Important: This feature will be only available from VASP 6.0.

Task

Calculation of the XANES K-edge in diamond using the supercell core-hole method.

Input

POSCAR

cubic diamond
 3.567
0.5 0.5 0.0
0.0 0.5 0.5
0.5 0.0 0.5
 2
direct
0.0 0.0 0.0
0.25 0.25 0.25

Above, we show the POSCAR file for the primitive unit cell of diamond. Note that we will not use this structure as input for the calculation. Instead, we use it to construct a POSCAR file for a super cell actually used in the calculation.

INCAR

System = DIAMOND
ALGO = FAST
ISMEAR = 0; SIGMA = 0.1;
ICORELEVEL = 2
CLNT = 1
CLN = 1
CLL = 0
CLZ = 1.0
CH_LSPEC = .TRUE.
CH_SIGMA = 0.5
NBANDS = 300
LREAL = A
  • To promote a core electron into the conduction bands and hence create the core-hole ICORELEVEL=2 has to be set. This corresponds to the final state approximation.
  • CLNT=1 selects the first atom species in the POSCAR file.
  • CLN=1 selects main quantum number 1 (hence K-edge).
  • CLL=0 selects angular quantum number 0 (s).
  • CLZ=1.0 selects the charge of the core hole. By setting this number to a fractional value we can mimic different screenings of the electrons. These non-integer values should only be used with caution, since this purely exploits error cancellation and does not correspond to a physically correct description of the screening.
  • By setting CH_LSPEC=.TRUE., we enable the calculation of matrix elements between core and conduction states and the calculation of the core electron absorption spectrum.
  • The broadening of the core electron absorption spectrum is controlled by the tag CH_SIGMA. Usually it is good practice to set this value low and broaden the spectrum in post processing.
  • We have to set NBANDS to a larger value to consider enough conduction band states in the calculation.
  • Since super cells are used the calculation of the projection operators in real space (LREAL=A) is much faster.

Calculation

Step 1 build a supercell

In the periodic boundary conditions of VASP, the core-hole interacts with its periodic replica so that we need sufficiently large super cells to reduce this spurious interaction. To this end, we employ successively large cells until the spectrum shows no significant changes. For this tutorial, we illustrate the calculation of a core-hole using a cell to allow for a reasonably fast calculation. However, for converged values one should use at least super cell.

The super cell can be obtained by either taking the file POSCAR.3x3x3 provided with this tutorial or constructing the POSCAR file from the primitive cell using p4vasp, which is demonstrated below:

  • Open p4vasp by typing p4v on the terminal.
  • Load the primitive cell by clicking on File→Load system:

  • Multiply cell in each direction (enter 3 for each direction) by clicking on Edit→Multiply Cell:

  • Save new system by clicking on File→Save system as:

Step 2 Prepare input files

The first few lines of the generated POSCAR file for the super cell should look like the following 

cubic diamond
3.567
 +1.5000000000  +1.5000000000  +0.0000000000
 +0.0000000000  +1.5000000000  +1.5000000000
 +1.5000000000  +0.0000000000  +1.5000000000
 54
Cartesian
 +0.0000000000  +0.0000000000  +0.0000000000
 +0.2500000000  +0.2500000000  +0.2500000000
 +0.5000000000  +0.0000000000  +0.5000000000
 ...

Here, all the 54 atoms are of the same species (line 6). To distinguish between the atom with the core-hole and the rest, we treat one atom as a different species. Choosing the first atom, we replace the 54 in the 6th line with 1 and 53 and obtain the following POSCAR file 

cubic diamond
3.567
 +1.5000000000  +1.5000000000  +0.0000000000
 +0.0000000000  +1.5000000000  +1.5000000000
 +1.5000000000  +0.0000000000  +1.5000000000
 1 53
Cartesian
 +0.0000000000  +0.0000000000  +0.0000000000
 +0.2500000000  +0.2500000000  +0.2500000000
 +0.5000000000  +0.0000000000  +0.5000000000
 ...

In the INCAR file, we specify that the first species carries the core-hole by setting CLNT=1. We create a POTCAR file with the PAW/PS information for both species. Since both species are carbon this amounts simply to the concatenation of the POTCAR file for carbon 

 cat POT_C POT_C > POTCAR

We provide the resulting POTCAR file for this core-hole calculation in the tar file of this tutorial.

To calculate accurate spectra, we need to include a sufficient number of conduction states. The required number of bands depends on the number of electrons in the system and the energy range of the spectrum. We can manually adjust the number of bands with the NBANDS variable. However, since the computation time increases drastically with the number of bands, selecting initially very large numbers is also not advisable. Hence, one has to increase the number of bands to find the optimum of computational effort and accuracy. In this tutorial, we use NBANDS=300.

The other input variables in the INCAR file are described above.

Mind: The multiplicity of the species carrying the core-hole has to be 1 otherwise the code will not work. Also mind that the selected species (CLNT in the INCAR file) is consistent with the order of the species specified in the POSCAR and POTCAR files.


Step 3 Running Calculation

Both the SCF calculation with the core-hole and the subsequent calculation of the dielectric matrix (spectrum) are done in a run of VASP. To minimize the spurious interaction between core-holes in neighboring cells requires large super cells and to reduce the computational time it is advisable to run a parallel VASP calculation 

 mpirun -np $np vasp_version

Here, $np corresponds to the number of processes and the _version in the executable stands for std, gam, and ncl: the standard, -point only, and non-collinear version, respectively. The cell used in this tutorial gives qualitatively correct results and can be completed even with a small number of processes. You can verify the spectrum on the larger cell, which we provide in the tar file of this tutorial.

Step 4 Extraction of XAS Spectrum

The XAS spectrum is proportional to the imaginary part of the frequency-dependent dielectric function, which is written in the OUTCAR file 

  frequency dependent IMAGINARY DIELECTRIC FUNCTION (independent particle, no local field effects) density-density
     E(ev)      X         Y         Z        XY        YZ        ZX
  --------------------------------------------------------------------------------------------------------------
  243.589609    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000
  243.677325    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000
  243.765042    0.000000    0.000000    0.000000    0.000000    0.000000    0.000000
  ...

Usually we are interested in the sum of all components of the dielectric matrix. You can obtain this by the script provided in this tutorial plot_core_imdiel.sh

Click to show/plot_core_imdiel.sh
#!/bin/bash

parallel=-1
normal=-1
all=-1
tauc=-1
trace=-1
while [[ $# -gt 0 ]]
do
   key="$1"
   case $key in
      -parallel) parallel=0
      ;;
      -normal) normal=0
      ;;
      -trace) trace=0
      ;;
      -tauc) tauc=0
      ;;
   esac
   shift
done


cat > helpscript.perl  <<EOF
#!/bin/perl

use strict;
use warnings;
my \$mode=shift;

while(<>)
{
   chomp;
   if(\$_ =~ /frequency dependent IMAGINARY DIELECTRIC FUNCTION/)
   {
      \$_=<>;
      \$_=<>;
      while (<>)
      {
         my \$sum=0;
         if (\$_ !~ /[0-9]/) {last;}
         chomp;
         \$_=~s/^/ /;
         my @help=split(/[\t,\s]+/);
         if (\$help[2]=~/NaN/||\$help[3]=~/NaN/||\$help[4]=~/NaN/) {next;}
         if (\$help[5]=~/NaN/||\$help[6]=~/NaN/||\$help[4]=~/NaN/) {next;}
         if (\$mode==0) {\$sum=\$help[2]+\$help[3]+\$help[4]+\$help[5]+\$help[6]+\$help[7];}
         if (\$mode==1) {\$sum=\$help[4];}
         if (\$mode==2) {\$sum=\$help[2]+\$help[3];}
         if (\$mode==3) {\$sum=\$help[2]+\$help[3]+\$help[4];}
         if (\$mode==4) {\$sum=(\$help[1]*\$help[1]*(\$help[2]+\$help[3]+\$help[4]+\$help[5]+\$help[6]+\$help[7]))**0.5;}
         if (\$mode==5) {\$sum=(\$help[1]*\$help[1]*(\$help[4]))**0.5;}
         if (\$mode==6) {\$sum=(\$help[1]*\$help[1]*(\$help[2]+\$help[3]))**0.5;}
         if (\$mode==7) {\$sum=(\$help[1]*\$help[1]*(\$help[2]+\$help[3]+\$help[4]))**0.5;}
         print \$help[1]," ",\$sum,"\n";
      }
   }
   last if eof;
}
EOF

if [[ $normal -eq 0 ]]; then
   if [[ $tauc -eq 0 ]]; then
      perl helpscript.perl 4 OUTCAR > CORE_DIELECTRIC_IMAG.dat
   else
      perl helpscript.perl 1 OUTCAR > CORE_DIELECTRIC_IMAG.dat
   fi
else
   if [[ $parallel -eq 0 ]]; then
      if [[ $tauc -eq 0 ]]; then
         perl helpscript.perl 5 OUTCAR > CORE_DIELECTRIC_IMAG.dat
      else
         perl helpscript.perl 2 OUTCAR > CORE_DIELECTRIC_IMAG.dat
      fi
   else
      if [[ $trace -eq 0 ]]; then
         if [[ $tauc -eq 0 ]]; then
            perl helpscript.perl 6 OUTCAR > CORE_DIELECTRIC_IMAG.dat
         else
            perl helpscript.perl 3 OUTCAR > CORE_DIELECTRIC_IMAG.dat
         fi
      else
         if [[ $tauc -eq 0 ]]; then
            perl helpscript.perl 7 OUTCAR > CORE_DIELECTRIC_IMAG.dat
         else
            perl helpscript.perl 0 OUTCAR > CORE_DIELECTRIC_IMAG.dat
         fi
      fi
   fi
fi
rm helpscript.perl

To use it type: 

bash ./plot_core_imdiel.sh

This will create the file CORE_DIELECTRIC_IMAG.dat containing the sum of the imaginary part of the dielectric matrix. Note that the absolute values of the experimental peak positions is not captured due to fundamental limitations of local DFT to describe the core level energies. Usually, there is even a noticeable deviation between calculations using different codes. Therefore, it is accepted in literature to look at the relative peak positions in the spectra and their relative intensity.

We compare the results obtained with VASP to experimental[1] and theoretical[2] XAS spectra from literature provided in the files C_XAS_aligned_to_VASP.dat and C_PARATEC_aligned_to_VASP.dat, respectively. The theoretical reference calculation was obtained using the PARATEC code and relies on PAW/Pseudopotential similar to VASP. We provide a script with this tutorial to compare these literature results to the spectrum obtained with VASP:

Click to show/gnuplot_XANES_C.script
unset ytics
set xrange [280:310]
set xlabel "Energy (eV)"
set ylabel "Absorption (arbitrary units)"
plot "CORE_DIELECTRIC_IMAG.dat" using 1:2 with lines lw 2 ti "VASP",\
     "C_PARATEC_aligned_to_VASP.dat" using 1:2 with lines ti "PAW lit",\
     "C_XAS_aligned_to_VASP.dat" using 1:2 with lines ti "Exp"
 pause -1

To use that script type: 

 gnuplot gnuplot_XANES_C.script

The file plot.sh constitutes a convenient wrapper around these post processing steps. The resulting spectra should look like this:

Because DFT cannot reproduce the absolute position (see above), we have shifted both spectra so that the position of the first peak coincides with the VASP. In this example, we scaled the experiment to VASP, since in this way the obtained results can be very easily compared using a script. Usually one would either scale the first peak to 0 or would scale the calculated value to the experiment. Additionally the intensity of the spectrum can be scaled arbitrarily. So in this example, we align the position and the height of the first peak for the calculations. The experiment is a little bit more tricky. It's a matter of taste what to consider as the first peak, but we decided that most likely the second peak corresponds to the first peak in the calculations and the first peak in experiment is a shoulder that is simply not pronounced in the calculations.

Another important issue is the broadening. Because the observed broadening is driven by many factors depending on the particular experimental setup, it is not possible to reproduce the broadening exactly. Therefore, we arbitrarily choose a broadening that gives approximately the same width as experiment. For simplicity in this tutorial, we choose to use a 0.5 eV constant Gaussian broadening (the broadening is determined by the ISMEAR tag and we used ISMEAR=0 which corresponds to a Gaussian broadening). For more elaborate spectra we strongly advise users to choose a 0.05 eV broadening and apply the desired broadening as post-processing.

Apart from the lower broadening width, we get a quite reasonable agreement with the theoretical literature calculation. We stress again that the super cell in this example are not fully converged. The interested user can repeat the calculations for a larger super cell. The files for this example are also given in the tar file. Be aware that the larger number of atoms in the bigger super cell requires an adjustment of the NBANDS variable..

Download

XANES_Diamond.tgz

References

  1. ↑ Y.Ma et al., Phys. Rev. Lett 69, 2598 (1992).
  2. ↑ M.Tallefumier et al., Phys. Rev. B 66, 195107 (2002).
Retrieved from "https://www.vasp.at/wiki/index.php?title=XANES_in_Diamond&oldid=17498"