# DFT-1/2 and DFT-PPS density functional methods for electronic structure calculations¶

**Version:** 2017.0

The 2017 release of QuantumATK introduces two novel density functional corrections for computing the electronic structure semiconductors and insulators: the DFT-1/2 and pseudopotential projector shift (DFT-PPS) methods. This tutorial gives an introduction to how to use these methods.

**DFT-1/2** is often also denoted LDA-1/2 or GGA-1/2,
and is a semi-empirical approach to correcting the self-interaction error
in local and semi-local exchange-correlation density functionals for extended systems,
similar in spirit to the DFT+U method.
As such, it can be viewed as an alternative to the TB09-Meta-GGA
method with broadly the same aims;
to improve the description of conduction-band energy levels and band gaps.
Just like with TB09-MGGA, the DFT+1/2 method is **not** suitable for calculations that
rely on total energy and forces, e.g., geometry optimization.

**PPS** is an abbreviation for pseudopotential projector shift.
This method introduces empirical shifts of the nonlocal projectors in
the SG15 pseudopotentials, in spirit of the empirical pseudopotentials proposed
by Zunger and co-workers in [WZ95]. The projector shifts have been
adjusted to reproduce technologically important properties of semiconductors
such as the fundamental band gap and lattice constant. The PPS method currently
works out-of-the-box for the elements silicon and germanium,
but may be used for other elements also by manually specifying the appropriate
projector shifts (these must somehow first be determined).
Importantly, the PPS method **does work** for geometry optimization
of semiconductor material structures.

## DFT-1/2 methods¶

The DFT-1/2 method attempts to correct the DFT self-interaction error by defining an atomic self-energy potential that cancels the electron-hole self-interaction energy. This potential is calculated for atomic sites in the system, and is defined as the difference between the potential of the neutral atom and that of a charged ion resulting from the removal of a fraction of its charge, between 0 and 1 electrons. The total self-energy potential is the sum of these atomic potentials. The addition of the DFT-1/2 self-energy potential to the DFT Hamiltonian has been found to greatly improve band gaps for a wide range of semiconducting and insulating systems [FMT11]. For more information, see the QuantumATK Manual entry on the DFT-1/2 method.

Important

The method is not entirely free of empirical parameters. For eaxmple, a cutoff radius is employed. Following Refs. [FMT08] and [FMT11] it is chosen such as to maximize the band gap.

Note also that not all elements in the system necessarily require the DFT-1/2 correction; it is generally advisable only to add this to the anionic species, and leave the cationic species as normal.

Default DFT-1/2 parameters are available in QuantumATK; these have been optimized against a wide range of materials, and should improve upon the standard LDA or GGA band gap in most cases.

### InP bandstructure using PBE-1/2¶

Create a new QuantumATK project (preferably in a new folder),
and open the **Builder**.
Then click
and locate the InP crystal in the **Database**.

Add the InP cystal to the Builder **Stash**, and send the configuration
to the **Script Generator**.
Then add the following script blocks:

and open the calculator widget to adjust the calculator parameters.

Select a 9x9x9 k-point grid, and navigate to the *Basis set/exchange
correlation* tab. Note that the default exchange-correlation method is
“GGA” (PBE). Enable the DFT-1/2 correction by ticking the checkbox
indicated below.

Close the calculator widget and save the script as `InP.py`

.
If you open the script using the **Editor**,
you will see that the exchange-correlation method is `GAHalf.PBE`

:

```
#----------------------------------------
# Exchange-Correlation
#----------------------------------------
exchange_correlation = GGAHalf.PBE
```

that is, GGA-1/2 using the PBE functional, which might also be denoted PBE-1/2.
You may download the script here if needed: `InP.py`

Run the script using either the **Job Manager**:

or from a terminal:

```
$ atkpython InP.py > InP.log
```

The calculation will be very fast.
Once done, use the **Bandstructure Analyzer** to plot the resulting band structure.
You will correctly find that InP is a direct-gap semiconductor
with a band band gap of 1.46 eV,
which is close to the experimental gap of 1.34 eV from Ref. [LRS96].

Tip

You may want to **zoom in** on the bands around the Fermi level for a better view.
There are two ways to do this:

- use the Zoom Tool (third icon from the left);
- or click the y-axis label to select it, then right-click and
select
*Edit item*in order to open the*Axis Properties*widget, where you can adjust the limits of the y-axis.

### III-V type semiconductor band gaps¶

The plots below show how the DFT-1/2 methods (LDA-1/2 and GGA-1/2) compare to LDA, GGA, and TB09-MGGA standard band gap calculations using default pseudopotentials and basis sets. A 9x9x9 k-point grid was used in all calculations, and experimental band gaps were adapted from Ref. [LRS96].

It is quite clear that the DFT-1/2 correction improves on standard
LDA and GGA band gaps, which are usually too small or non-existing
(no bar). TB09-MGGA calculations with self-consistently determined *c*-parameter
also improves band gaps in general.

The red bars for Ge using LDA-1/2 and GGA-1/2 indicate that a *direct*
band gap is wrongly predicted (\(\Gamma\)–\(\Gamma\) instead of \(\Gamma\)–L
as in experiments).
The orange bar for GaP using GGA-1/2 indicates that the gap is correctly
predicted to be *indirect*, but that the band energy in the L-valley
is wrongly predicted to be a little lower than in the X-valley.

### Manually specifying DFT-1/2 parameters¶

It is of course possible to manually specify the DFT-1/2 parameters
instead of using the default ones.
This is done by creating an instance of the `DFTHalfParameters`

class, which is then given as an argument to the basis set.

For example, let’s consider the case of GaAs. The following script manually sets up
the As DFT-1/2 parameters identical to the default ones,
and leaves the DFT-1/2 correction disabled for Ga, which is also default
behavior. Note the `dft_half_parameters`

argument to the basis set:

```
#----------------------------------------
# Basis Set
#----------------------------------------
# LDA-1/2 parameters for As
dft_half_parameters = DFTHalfParameters(
element=Arsenic,
fractional_charge=[0.3, 0.0],
cutoff_radius=4.0*Bohr,
)
# No LDA-1/2 parameters are needed for Ga (Disabled)
basis_set = [
LDABasis.Arsenic_DoubleZetaPolarized(
dft_half_parameters=dft_half_parameters),
LDABasis.Gallium_DoubleZetaPolarized(
dft_half_parameters=Disabled),
]
```

For more information, see the sections Input format for DFT-1/2 and DFT-1/2 default parameters in the QuantumATK Manual.

Warning

Choosing the appropriate DFT-1/2 parameters may be a very delicate matter, and great care has been taken in determining the default parameters. If you choose to use non-default DFT-1/2 parameters, the quality of those parameters is entirely your own responsibility as a user!

QuantumWise does not offer support for determining custom DFT-1/2 parameters; we recommend in general that users stick to the default ones.

## DFT-PPS method¶

As already mentioned, the DFT-PPS method applies shifts to the nonlocal projectors in the SG15 pseudopotentials. The nonlocal part of the pseudopotential, \(\hat{V}_\text{nl}\), is modified according to

where the sum is over all projectors \(p_{l}\), and \(\alpha_{l}\) is an empirical parameter that depends on the orbital angular momentum quantum number \(l\). Note that this approach does not increase the computational cost of DFT calculations!

The required projector shift parameters have been optimized for silicon and germanium and for use with the PBE density functional and SG15 pseudopotentials only. These are implemented as separate basis sets:

```
BasisGGASG15.Silicon_LowProjectorShift
BasisGGASG15.Silicon_MediumProjectorShift
BasisGGASG15.Silicon_HighProjectorShift
BasisGGASG15.Silicon_UltraProjectorShift
BasisGGASG15.Germanium_MediumProjectorShift
BasisGGASG15.Germanium_HighProjectorShift
BasisGGASG15.Germanium_UltraProjectorShift
```

For each element, the same set of optimized projector shifts are applied in
all SG15 basis sets. The script `projector_shifts.py`

prints the built-in DFT-PPS parameters
for Si and Ge:

```
basis_sets = [
BasisGGASG15.Silicon_MediumProjectorShift, # Si PPS-PBE SG15-Medium
BasisGGASG15.Germanium_HighProjectorShift , # Ge PPS-PBE SG15-High
]
for basis_set,element in zip(basis_sets,['Si','Ge']):
print(element)
projector_shift = basis_set.projectorShift()
print("s-shift: %+.3f eV" % projector_shift.sOrbitalShift().inUnitsOf(eV))
print("p-shift: %+.3f eV" % projector_shift.pOrbitalShift().inUnitsOf(eV))
print("d-shift: %+.3f eV" % projector_shift.dOrbitalShift().inUnitsOf(eV))
```

Running it produces the output shown below. The `d-shift`

for Si is 0.0 eV
since silicon has no *d*-electrons:

```
Si
s-shift: +21.330 eV
p-shift: -1.430 eV
d-shift: +0.000 eV
Ge
s-shift: +13.790 eV
p-shift: +0.220 eV
d-shift: -2.030 eV
```

### Si, SiGe, and Ge band gaps and lattice constants¶

The DFT-PPS method is enabled in the **Script Generator**
by selecting one of the `ProjectorShift`

basis sets for Si or Ge,
as indicated below. No other calculator settings need to be changed in order
to switch from ordinary PBE to DFT-PPS.

One very convenient aspect of the DFT-PPS method is that it allows for geometry optimization (forces and stress minimization) just like ordinary GGA calculations do – in fact, the DFT-PPS parameters may often be chosen such as to give highly accurate semiconductor lattice constants while also producing accurate band gaps.

In the following, you will consider bulk Si and Ge, and a simple 50/50
SiGe alloy. These 3 bulk configurations are defined in the script `bulks.py`

.
The script `pbe.py`

runs geometry optimization and band structure analysis
for all 3 configurations,
while the script `pps.py`

does the same using the DFT-PPS method.

The latter script looks like shown below. Note the lines from 2 to 12 in that particular script, where the bulk configurations are imported from the external script, and a Python loop over these configurations and basis sets is set up:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 | ```
# -*- coding: utf-8 -*-
from bulks import si, ge, sige
setVerbosity(MinimalLog)
configurations = [si,ge,sige]
labels = ['Si','Ge','SiGe']
si_basis = BasisGGASG15.Silicon_MediumProjectorShift
ge_basis = BasisGGASG15.Germanium_HighProjectorShift
basis_sets = [[si_basis], [ge_basis], [si_basis,ge_basis]]
for bulk_configuration,label,basis_set in zip(configurations,labels,basis_sets):
outfile = "%s_PPS.hdf5" % label
# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------
k_point_sampling = MonkhorstPackGrid(
na=9,
nb=9,
nc=9,
)
numerical_accuracy_parameters = NumericalAccuracyParameters(
k_point_sampling=k_point_sampling,
density_mesh_cutoff=90.0*Hartree,
)
calculator = LCAOCalculator(
basis_set=basis_set,
numerical_accuracy_parameters=numerical_accuracy_parameters,
)
bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()
nlsave(outfile, bulk_configuration)
# -------------------------------------------------------------
# Optimize Geometry
# -------------------------------------------------------------
fix_atom_indices_0 = [0, 1]
constraints = [FixAtomConstraints(fix_atom_indices_0)]
bulk_configuration = OptimizeGeometry(
bulk_configuration,
max_forces=0.05*eV/Ang,
max_stress=0.1*GPa,
max_steps=200,
max_step_length=0.2*Ang,
constraints=constraints,
trajectory_filename=None,
optimizer_method=LBFGS(),
constrain_bravais_lattice=True,
)
nlsave(outfile, bulk_configuration)
nlprint(bulk_configuration)
# -------------------------------------------------------------
# Bandstructure
# -------------------------------------------------------------
bandstructure = Bandstructure(
configuration=bulk_configuration,
route=['L', 'G', 'X'],
points_per_segment=50
)
nlsave(outfile, bandstructure)
``` |

Download all 3 scripts (`bulks.py`

, `pbe.py`

, `pps.py`

) and run the latter two
using either the **Job Manager** or from a terminal:

```
$ atkpython pps.py > pps.log
$ atkpython pbe.py > pbe.log
```

The jobs should take around 5 minutes each.

Use then the script `plot_pps.py`

to plot the results.
Running that script should produce the figure shown below,
where the calculated indirect band gaps are shown in red circles,
and the calculated lattice constants are shown in blue squares,
both for PBE (dashed lines) and DFT-PPS (solid lines).

The DFT-PPS band gaps of Si and Ge compare very well to experiments (block dots; from Ref. [LRS96]), and varies roughly linearly with germanium content. On the contrary, the ordinary PBE method predicts zero Ge band gap.

The DFT-PPS lattice constants of pure Si and Ge are also closer to experiments (grey squares) than found with non-corrected PBE.

### Manually specifying DFT-PPS parameters¶

It is of course possible to manually specify the DFT-PPS projector shift parameters instead of using the default ones. This may be especially useful in DFT-PPS calculations for elements where there are no default DFT-PPS parameters (only Si and Ge currently have defaults).

The projector shifts are supplied using an instance of the
`PseudoPotentialProjectorShift`

class,
which is then given as an argument to the SG15 basis set.
As an example, the following script manually sets the Si and Ge DFT-PPS
parameters identical to the default ones:

```
#----------------------------------------
# Basis Set
#----------------------------------------
# Basis set for Silicon
SiliconBasis_projector_shift = PseudoPotentialProjectorShift(
s_orbital_shift=21.33*eV,
p_orbital_shift=-1.43*eV,
d_orbital_shift=0.0*eV,
f_orbital_shift=0.0*eV,
g_orbital_shift=0.0*eV
)
SiliconBasis = BasisGGASG15.Silicon_Medium(projector_shift=SiliconBasis_projector_shift)
# Basis set for Germanium
GermaniumBasis_projector_shift = PseudoPotentialProjectorShift(
s_orbital_shift=13.79*eV,
p_orbital_shift=0.22*eV,
d_orbital_shift=-2.03*eV,
f_orbital_shift=0.0*eV,
g_orbital_shift=0.0*eV
)
GermaniumBasis = BasisGGASG15.Germanium_High(projector_shift=GermaniumBasis_projector_shift)
# Total basis set
basis_set = [
SiliconBasis,
GermaniumBasis,
]
```

Warning

Choosing the appropriate DFT-PPS parameters may be a very delicate matter, and often requires a numerical optimization procedure. The SciPy software offers numerous such routines, but if you choose to use non-default DFT-PPS parameters, the quality of those parameters is entirely your own responsibility as a user!

QuantumWise does not offer support for optimizing DFT-PPS parameters; we recommend in general that users stick to the default parameters, if they exist.

## References¶

[FMT08] | Luiz G. Ferreira, Marcelo Marques, and Lara K. Teles. Approximation to density functional theory for the calculation of band gaps of semiconductors. Phys. Rev. B, 78:125116, Sep 2008. doi:10.1103/PhysRevB.78.125116. |

[FMT11] | (1, 2) Luiz G. Ferreira, Marcelo Marques, and Lara K. Teles. Slater half-occupation technique revisited: the LDA-1/2 and GGA-1/2 approaches for atomic ionization energies and band gaps in semiconductors. AIP Adv., 1(3):032119, 2011. doi:10.1063/1.3624562. |

[LRS96] | (1, 2, 3) M. Levinshtein, S. Rumyantsev, and M. Shur, editors. Handbook Series on Semiconductor Parameters. volume 1. World Scientific Publishing Cp. Pte. Ltd., Singapore, 1996. |

[WZ95] | L.-W. Wang and A. Zunger. Local-density-derived semiempirical pseudopotentials. Phys. Rev. B, 51:17398–17416, 1995. doi:10.1103/PhysRevB.51.17398. |