PBE and HSE06 band structures of GaAs

Version: 2016.4

In this tutorial you will learn the basics of how to use the FHI-aims calculator in ATK. FHI-aims is an all-electron electronic structure simulations package, so it does not employ any pseudopotentials. Moreover, it provides access to a wide range of beyond-GGA methods, including hybrid density functionals, the random-phase approximation for total energies, and the GW approximation for excited states. See the FHI-aims website for more details.

Attention

Due to the Synopsys acquisition of QuantumWise in the fall of 2017, see the Synopsys press release, it is no longer possible to obtain a VNL-ATK license that includes access to the FHI-aims calculator.

introbar

Introduction

We shall consider the band structure of GaAs using the PBE and HSE06 density functionals, the latter of which employs a screened variant of exact exchange. The HSE06 can be computationally intense, but is often more accurate than standard GGA methods. See also Refs. [HSE03], [PMH+06], and [KVIS06].

Tip

The video below considers silicon as an example semiconductor, while the written below it gives instructions using GaAs instead. Do not let this confuse you; the difference in workflow is very small.

ATK 2017: The video was made before the default ATK datafile format was changed to HDF5. Therefore, ATK 2017 users should just use the ”.hdf5” extension wherever the tutorial text says ”.nc”.

Many more video tutorials are available on the QuantumWise YouTube channel.

Note

Your ATK license file must include the ATKFHIaims feature in order to use FHI-aims with ATK. Please contact us for any inquiries.

Important

This tutorial uses rather low accuracy settings in order to keep the computational time to a minimum, for instructional purposes. Please remember that settings related to k-point grids, basis sets, and integration grids should always be carefully checked for convergence in publication-quality scientific work.

PBE band structure

  • Open the builder_icon Builder and set up the GaAs crystal (zincblende structure).
  • Send the structure to the script_generator_icon Scripter and add the calculator_icon New Calculator and analysis_icon Analysis ‣ Bandstructure script blocks.
  • Open the New Calculator block and select the FHI-aims calculator.
  • Select the PBE density functional and set up a 12x12x12 k-point grid.
  • In Control Parameters, use the system default settings for “Insulator”.
../../_images/pic_20.png

Save the script as pbe.py (but do not close the Script Generator window) and run the script using the job_manager_icon Job Manager. Alternatively, you can run the FHI-aims calculation from the command line:

atkpython pbe.py > pbe.log

The calculation should take only a few minutes, and the GaAs band structure will pop up on the VNL LabFloor. You can visualize it using the Bandstructure Analyzer and zoom in around the Fermi level:

../../_images/pic_21.png

HSE06 band structure

You will now compute the same band structure using HSE06 and compare it to the PBE results. HSE06 can be computationally heavy for periodic systems, so we will here use a coarser 6x6x6 k-point grid.

Note

In the following, the PBE and HSE06 calculations are conveniently done after each other in a single Python script. This illustrates one of the advantages of running FHI-aims calculations through ATK Python scripts: Composite workflows are easily created using the flexibility of Python scripting, see more at Python in ATK.

  • PBE calculation:
    • Go back to the script_generator_icon Script Generator window used above. Set the k-point grid to 6x6x6, and change the band structure route to L, G, X.
  • HSE calculation:
    • Add one more pair of calculator_icon New Calculator and analysis_icon Bandstructure blocks.
    • Adjust the settings to make them identical to those for PBE, but chooce HSE06 as xc functional.

Finally, change the default output file name to hse.nc, and save the script as hse.py.

../../_images/hse.png

Run the calculations, preferably in parallel using e.g. 4 MPI processes:

../../_images/jm.png

You can of course obtain the same in a terminal:

QWDIR/libexec/mpiexec.hydra -n 4 atkpython hse.py > hse.log

where QWDIR is the VNL-ATK installation directory.

The calculations should take roughly 30 minutes in total if parallelized over 4 CPU cores on a modern laptop.

When finished, you can mark both band structure objects on the LabFloor and use the Compare Data plugin to visualize both band structures in one plot:

../../_images/pic_31.png ../../_images/pic_30.png

The two band structures are clearly different; the HSE06 conduction bands (green, object ID gID003) are pushed up relative to the PBE ones, leading to a predicted band gap of 1.36 eV, not too far from the experimental value of 1.52 eV.

Note

Please be aware that the results above were obtained using Medium accuracy (“Light” integration grids and “Tier 1” basis sets) and a somewhat coarse k-point grid. Calculations with Tight settings, more basis functions, and more k-points may give different results. See more below.

Also note that band structure calculations should usually be preceeded by a first principles geometry optimization of the structure.

Basis sets and integration grids

../../_images/GaAs.png

Fig. 105 GaAs band gap calculated with different FHI-aims settings.

As seen in the figure above, FHI-aims results do in general depend on the computational settings, for example basis sets and integration grids. “Tight” means more refined integrations grids than “Light”, while “Tier 2” has more basis functions than “Tier 1”. Note that increased computational quality requires in general also more computing resources (CPU hours and/or memory).

The calculations were performed using the script GaAs.py, which is reproduced below. The first 18 script lines (highlighted) set up a Python loop over different combinations of basis sets and levels of integration accuracy. SCF and band structure calculations for all 6 combinations are hereby executed in a single ATK job.

# ====================================================
basis_sets  = [Light, Tight]
tiers       = [0, 1, 2]
# ====================================================

for basis in basis_sets:
    for tier in tiers:
        out = 'GaAs_pbe_%s_tier-%i.nc' % (basis.speciesDefaultAsString(), tier)

        #----------------------------------------
        # Basis Set
        #----------------------------------------
        basis_set = FHIaimsBasisSet()
        species_bases = [
            FHIaimsBasisSetSpecies(Arsenic, basis, tier),
            FHIaimsBasisSetSpecies(Gallium, basis, tier),
            ]
        basis_set.setAndCheckSpeciesBasis(species_bases)

        # -------------------------------------------------------------
        # Bulk Configuration
        # -------------------------------------------------------------
        # Set up lattice
        lattice = FaceCenteredCubic(5.6537*Angstrom)
        # Define elements
        elements = [Gallium, Arsenic]
        # Define coordinates
        fractional_coordinates = [[ 0.  ,  0.  ,  0.  ],
                                  [ 0.25,  0.25,  0.25]]
        # Set up configuration
        bulk_configuration = BulkConfiguration(
            bravais_lattice=lattice,
            elements=elements,
            fractional_coordinates=fractional_coordinates
            )

        # -------------------------------------------------------------
        # Calculator
        # -------------------------------------------------------------
        k_point_sampling = MonkhorstPackGrid(9,9,9)
        numerical_accuracy_parameters = NumericalAccuracyParameters(
            k_point_sampling=k_point_sampling,
            )
        iteration_control_parameters = IterationControlParameters(
            damping_factor=0.6,
            algorithm=PulayMixer(),
            )

        #----------------------------------------
        # FHI-aims Parameters
        #----------------------------------------
        fhiaims_parameters = FHIaimsParameters()
        fhiaims_parameters.addControlKeywordValuePair("sc_accuracy_eev","0.001")
        fhiaims_parameters.addControlKeywordValuePair("charge","0")
        fhiaims_parameters.addControlKeywordValuePair("sc_accuracy_etot","1e-06")
        fhiaims_parameters.addControlKeywordValuePair("sc_accuracy_rho","1e-05")
        fhiaims_parameters.addControlKeywordValuePair("occupation_type","gaussian 0.000001")

        #----------------------------------------
        # Working Directory
        #----------------------------------------
        working_directory = '.'

        #----------------------------------------
        # FHI-aims calculator
        #----------------------------------------
        calculator = FHIaimsCalculator(
            basis_set=basis_set,
            exchange_correlation=FHIaims_GGA_PBE,
            numerical_accuracy_parameters=numerical_accuracy_parameters,
            iteration_control_parameters=iteration_control_parameters,
            fhiaims_parameters=fhiaims_parameters,
            working_directory=working_directory,
            )

        bulk_configuration.setCalculator(calculator)
        nlprint(bulk_configuration)
        bulk_configuration.update()

        # -------------------------------------------------------------
        # Bandstructure
        # -------------------------------------------------------------
        bandstructure = Bandstructure(
            configuration=bulk_configuration,
            route=['G','L'],
            points_per_segment=5,
            )
        nlsave(out, bandstructure)

The script plot.py may then be used to read the band structures, extract the band gap, plot the results, and save the plot as a PNG file:

# READING DATA
basis_sets  = [Light, Tight]
tiers       = [0, 1, 2]
gaps   = []
labels = []
for tier in tiers:
    for basis in basis_sets:
        out = 'GaAs_pbe_%s_tier-%i.nc' % (basis.speciesDefaultAsString(), tier)
        bandstructure = nlread(out, Bandstructure)[0]
        gap = bandstructure._directBandGap().inUnitsOf(eV)
        gaps.append(gap)
        label = "%s\n" % basis.speciesDefaultAsString().capitalize()
        label += "Tier %i" % tier
        labels.append(label)

# PLOTTING DATA
import pylab
x = range(len(labels))
pylab.figure(figsize=(7,4))
pylab.plot(x, gaps, 'bo-')
pylab.xlim((-0.2,5.2))
pylab.ylabel('GaAs direct band gap (eV)')
pylab.title('XC=PBE using 9x9x9 k-points')
ax = pylab.gca()
ax.set_xticks(x)
ax.set_xticklabels(labels)
pylab.tight_layout()
pylab.savefig('GaAs.png')
pylab.show()

You can try to run the scripts yourself – the calculations should take ca. 30 mins in total if run in serial on a modern laptop.

Germanium band structure

For completeness, you can also calculate the band structure of Ge using HSE06: Download the script germanium.py and run the calculation. It takes roughly 25 minutes when parallelized on 16 cores, but can of course be run on fewer cores. The result is illustrated below.

../../_images/hse_germanium.png

References

[HSE03]Jochen Heyd, Gustavo E. Scuseria, and Matthias Ernzerhof. Hybrid functionals based on a screened coulomb potential. J. Chem. Phys., 118(18):8207–8215, 2003. doi:10.1063/1.1564060.
[KVIS06]Aliaksandr V. Krukau, Oleg A. Vydrov, Artur F. Izmaylov, and Gustavo E. Scuseria. Influence of the exchange screening parameter on the performance of screened hybrid functionals. J. Chem. Phys., 2006. doi:10.1063/1.2404663.
[PMH+06]J. Paier, M. Marsman, K. Hummer, G. Kresse, I. C. Gerber, and J. G. Angyan. Screened hybrid density functionals applied to solids. J. Chem. Phys., 2006. doi:10.1063/1.2187006.