Hybrid-Functional Calculations¶

by Cecilia Vona, Ute Werner, & Dmitrii Nabok for exciting neon

(Jupyter notebook by Mara Voiculescu & Martin Kuban)

Purpose: In this tutorial, you will learn how to perform ground-state calculations with the hybrid functionals PBE0 and HSE. We will use diamond as test material. Since Diamond has an indirect band-gap, instead of considering the gap, we will look at how the transition energies in specific symmetry point changes, depending on the functional employed. Specifically we will consider the transition Γ→Γ and Γ→X.

Table of Contents

0. Before Starting

1. Theoretical Background

2. PBE Calculation

Preparation of the Input File
Execute the Calculation

3. PBE0 Calculation

Preparation of the Input File
Execute the Calculation

4. HSE Calculation

Preparation of the Input File
Execute the Calculation

5. Post Processing

6. Band Structure and Density of States

7. Exercise

0. Before Starting¶

Read the following paragraphs before starting with the rest of this tutorial!

Before running any Jupyter tutorials, please refer to the 00_before_starting.md document on how to correctly set up the environment. This only needs to be done once. After which, the venv can be (re)activated from exciting's root directory:

source tools/excitingjupyter/venv/excitingvenv/bin/activate

As a first step, you may create a running directory for the notebook.

%%bash
mkdir -p run_tutorial_hybrid_functional_calculations

1. Theoretical Background¶

The hybrid exchange-correlation functionals are constructed by substituting a fraction of local/semi-local exchange by a fraction of non-local exact-exchange:

\begin{equation} E^{hyb}_{xc} = E^L_{xc}+\alpha (E^{NL}_x-E^L_x). \end{equation}

The non-local exact exchange is often computed employing Hartree-Fock (HF). The type of the DFT local/semi-local functional and the fraction of non-local exact exchange substituted ($\alpha$) vary depending on the hybrid functional type.

The hybrid functionals PBE0 and HSE are implemented in exciting. PBE0 has the following form:

\begin{equation} E^{PBE0}_{xc} = E^{PBE}_{xc}+\alpha (E^{HF}_x-E^{PBE}_x)\; , \end{equation}

in which the semi-local functional is PBE and the fraction of exact-exchange substituted is uniquely determined by the mixing parameter α, which has the value 0.25.

The hybrid functional HSE, which is a screened hybrid functional, has the following form:

\begin{equation} E^{HSE}_{xc} = E^{PBE}_{xc}+\alpha [E^{HF,SR}_x(\omega)-E^{PBE,SR}_x(\omega)]. \end{equation}

Also in HSE, PBE is used as semi-local functional, but only the short-range part of the non-local exact exchange is substituted. In fact, in HSE the amount of non-local exact exchange is determined by both the mixing parameter $\alpha$=0.25 and the screening screening parameter $\omega$, which defines the cutoff between the short-range and the long-range part of the Coulomb operator through the error function and its complementary:

\begin{equation} v(r)=v^{SR}(r,\omega)+v^{LR}(r,\omega)=\frac{\text{erfc}(\omega r)}{r}+\frac{\text{erf}(\omega r)}{r}. \end{equation}

The version of HSE which is implemented in exciting is HSE06, in which the PBE and HF contributions have the same value of $\omega$=0.11 $a_0^{-1}$.

2. PBE Calculation¶

In order to investigate and compare the effect of hybrid functionals on the transition energies, we will start by performing a GGA calculation. For this purpose, we will use the PBE functional.

Before starting the calculation, create a subdirectory called PBE.

%%bash
cd run_tutorial_hybrid_functional_calculations
mkdir -p PBE
cd ..

i) Preparation of the Input File¶

We start by creating an exciting (xml) input file called input.xml to compute the ground-state of diamond using the semi-local functional PBE. The file can look like this:

<input>

   <title>Diamond PBE</title>

   <structure speciespath="$EXCITINGROOT/species" >
      <crystal scale="6.7425">
         <basevect> 0.0     0.5     0.5 </basevect>
         <basevect> 0.5     0.0     0.5 </basevect>
         <basevect> 0.5     0.5     0.0 </basevect>
      </crystal>
      <species speciesfile="C.xml">
         <atom coord="0.00 0.00 0.00" />
         <atom coord="0.25 0.25 0.25" />
      </species>
   </structure>

   <groundstate 
      ngridk="4 4 4"
      rgkmax="5.0"
      xctype="GGA_PBE">
   </groundstate>

</input>

The next step is writing the complete input as a string and saving it in your working directory as input.xml.

ii) Execute the Calculation¶

Before starting the calculation, remember to set the correct path for the exciting root directory ($EXCITINGROOT) to the one pointing to the place where the exciting directory is placed, using the command

%%bash
cd run_tutorial_hybrid_functional_calculations/PBE
python3 -m excitingscripts.setup.excitingroot
cd ../..

Start the calculation by using the script excitingscripts.execute.single.

%%bash
cd run_tutorial_hybrid_functional_calculations
python3 -m excitingscripts.execute.single -r PBE
cd ..

3. PBE0 Calculation¶

Now we can repeat the calculation by using hybrid functionals. We will start with the hybrid functional PBE0.

i) Preparation of the Input File¶

To compute the ground state of diamond with the hybrid functional PBE0 it is necessary to modify the groundstate element and, if you like, you can modify/adapt the content of the title element as in the example

<input>

<title>Diamond PBE0</title>

...   

   <groundstate 
      ngridk="4 4 4"
      rgkmax="5.0"
      nempty="20"
      xctype="HYB_PBE0">
      <Hybrid  
         excoeff="0.25" />
   </groundstate>

</input>

The first thing to notice is the presence of the new element Hybrid inside the element groundstate. This element is not mandatory to perform calculations with hybrid functionals. If it is not included, all the parameters specified by the element Hybrid will be initialized with their default values. Moreover, by comparing the input with the input of PBE, we observe that in the example is specified the number of empty states (nempty), which is crucial for having convergent hybrids calculations.

To the parameters introduced in the ground-state calculation with hybrid functionals will be discussed later. Before we do that, execute the calculation with PBE0 and create the HSE input file.

%%bash
cd run_tutorial_hybrid_functional_calculations
mkdir -p PBE0
cd ..

%%bash
cd run_tutorial_hybrid_functional_calculations/PBE0
python3 -m excitingscripts.setup.excitingroot
cd ../..

ii) Execute the Calculation¶

%%bash
cd run_tutorial_hybrid_functional_calculations
python3 -m excitingscripts.execute.single -r PBE0
cd ..

4. HSE Calculation¶

Now we repeat everything using the hybrid functional HSE.

i) Preparation of the Input File¶

Start by modifying again the groundstate element in the file input.xml, inside the directory diamond-hybrids, as in the example:

<input>

<title>Diamond HSE</title>

... 

   <groundstate 
      ngridk="4 4 4"
      rgkmax="5.0"
      nempty="20"
      xctype="HYB_HSE">
      <Hybrid 
         excoeff="0.25" 
         omega="0.11" />
   </groundstate>

</input>

%%bash cd run_tutorial_hybrid_functional_calculations mkdir -p HSE cd ..

%%bash
cd run_tutorial_hybrid_functional_calculations/HSE
python3 -m excitingscripts.setup.excitingroot
cd ../..

ii) Execute the Calculation¶

%%bash
cd run_tutorial_hybrid_functional_calculations
python3 -m excitingscripts.execute.single -r HSE
cd ..

You have probably noticed, that the running time slightly increases when HSE is employed instead of PBE0. This can seem weird, since HSE is known for being faster. For a fix number of k-points this is not the case in exciting, due to the LAPW+lo basis. However, convergent HSE calculations are usually faster than convergent PBE0 calculations, since the latter converge considerably slower with respect to the number of k-points.

$\Rightarrow$ Plot k-points convergence

These plots show the difference of the transition energies between Γ→Γ and Γ→X, computed with the PBE, PBE0 and HSE xc-functionals, for different k-points grid. The reference calculation is computed for $N_k=10^3$. From this you can have a better understanding of the k-points convergence of the different functionals.

To repeat the calculations with PBE0 and HSE, for higher number of k-points, requires too much time for the purpose of this tutorial. You can try to reproduce these plots as an exercise.

We can now have a closer look to the most relevant parameters in ground-state calculations with hybrid functionals.

Parameter	Description
`xctype`="HYB_PBE0" or "HYB_HSE"	The type of (hybrid) functional used needs to be specified in `xctype`. In `exciting` the hybrid functionals available are PBE0 or HSE and the one shown are the respectively keywords.
`nempty`="20"	The number of empty states is crucial for the convergence and is dependent from `rgkmax`.
`excoeff`="0.25"	Keyword, inside the element `Hybrid`, for the mixing parameter $\alpha$. The default value is 0.25.
`omega`="0.11"	Keyword, inside the element `Hybrid`, for to screening parameter $\omega$ in units of $a_0^{-1}$ (only with HSE). The default value is 0.11 $a_0^{-1}$.

For further details on the parameters see Input Reference.

5. Post Processing¶

To observe some of the effects that the different functionals produce on the electronic structure, we compare the transition energies in specific high-symmetry points. Specifically, we will consider the following transition Γ→Γ and Γ→X.

To do this, move into the parent directory and run the script excitingscripts.compare_transition_energies.

%%bash
cd run_tutorial_hybrid_functional_calculations
python3 -m excitingscripts.compare_transition_energies -r PBE PBE0 HSE
cd ..

The output should look like this:

------------------------------------------------

 Transition energies in eV:

                   PBE                PBE0               HSE                
Gamma -> Gamma:    5.631773175190849  7.908401189753921  7.121367501481966  
Gamma -> X:        4.811783049188941  6.915520695999625  6.120261918087781  

------------------------------------------------`

6. Band Structure and Density of States¶

When hybrids functionals are employed, to compute properties, such as band structure and density of states, the Wannier-functions interpolation schema needs to be adopted. To use the Wannier interpolation schema is not straightforward and therefore there is a tutorial dedicated Wannier Functions for Interpolation in Reciprocal Space.

In the tutorials the properties are computed only for the hybrid functional PBE0, as exercise you can try to repeat the Wannier-functions tutorial for the hybrid functional HSE.

7. Exercise¶

Try to modify some of the parameters such as rgkmax, nempty or ngridk in the input files. Remember that the optimal value of nempty depends on the value of rgkmax. For more complex systems nempty can be quite big.

In order to get an hint about the maximum number of empty states that can be considered, use the following command.

grep "Maximum Hamiltonian size" INFO.OUT

Literature¶

C. Adamo and V. Barone, J. Chem. Phys. 110, 6158, (1999)
M. Ernzerhof and G. E. Scuseria, J. Chem. Phys. 110, 5029 (1999)
A. V. Krukau et al., J. Chem. Phys. 125, 224106 (2006)
M. Schlipf et al., Phys. Rev. B 84, 125142 (2011)