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

Version: 2017.0

The 2017 release of ATK 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 Meta-GGA method with broadly the same aims; to improve the description of conduction-band energy levels and band gaps. Just like with 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 ATK 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 ATK; 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 VNL project (preferably in a new folder), and open the Builder. Then click Add ‣ From Database 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:

• New Calculator
• Bandstructure

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: 1. use the Zoom Tool (third icon from the left); 2. 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 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). 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 ATK 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 $\hat{V}_\text{nl} \rightarrow \hat{V}_\text{nl} \mathrel{+} \sum_l |p_{l} \rangle \alpha_{l} \langle p_{l} | ,$ 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.