Warning

Requires optional libraries

This notebook requires a working installation of GPAW.

conda

Install GPAW using conda:

conda install -c conda-forge gpaw

pip

Install GPAW using pip:

pip install gpaw

Install the PAW datasets into the folder <dir> using this command:

gpaw install-data <dir>

See here for detailed installation instructions.

Show code cell source Hide code cell source

import abtem
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from ase.build import graphene

from gpaw import GPAW

abtem.config.set({"local_diagnostics.progress_bar": False});

ab initio potentials with GPAW

The independent atom model (IAM) neglects bonding and the resulting charge transfer, this is typically a good approximation, because the core charge and electrons are responible for most of the potential. However, the difference between the IAM and a more realistic model including the effects of bonding is measurable, in certain systems[MPS21]. Understanding this charge transfer may provide further insights, due its the importance in chemistry.

Here we go beyond the independent atom model using density functional theory (DFT), to calculate a potential that includes charge transfer using GPAW[MHJ05]. Note that you need a working GPAW installation to run the code in this walkthrough, see the GPAW documentation for more information. abTEM only supports GPAW directly, however, other DFT codes may be used assuming that you can obtain the electron charge density as a numpy array through the ChargeDensityPotential.

DFT calculation of hBN with GPAW

The first step to create a DFT potential is to converge a DFT calculation. We calculate the minimal hexagonal cell of hexagonal-Boron-Nitride (hBN). The atomic model is created by modifying a graphene model.

atoms = graphene(vacuum=3, a=2.504)

atoms[0].number = 5
atoms[1].number = 7

abtem.show_atoms(atoms, legend=True);

To run the DFT calculation, we use the GPAW calculator. We use the default parameters, except for the Brillouin zone sampling, which for a cell this small should be at least \(5 \times 5 \times 3\).

gpaw = GPAW(txt=None, kpts=(5, 5, 1))
atoms.calc = gpaw

Running the method get_potential_energy triggers the DFT self-consistent field cycle to run, after which the GPAW object contains the converged PAW electron density.

atoms.get_potential_energy()

/home/jacob/miniconda3/envs/abtem/lib/python3.11/site-packages/gpaw/calculator.py:713: DeprecatedParameterWarning: Finite-difference mode implicitly chosen; it will be an error to not specify a mode in the future
  warnings.warn(

-19.36957650208433

Warning

DFT calculations can be extremely computationally intensive and may require massive parallelization. Hence, it may not be appropriate to run the DFT simulation in a notebook. Instead, we recommend exporting the converged GPAW calculator to a file, then importing it for your abTEM simulation.

Using the GPAW potential in abTEM

It is straightforward to calculate a DFT potential from a converged GPAW calculation. The GPAWPotential object just requires a converged GPAW calculator.

potential_dft = abtem.GPAWPotential(gpaw, sampling=0.04).build().compute()

The GPAWPotential shares the same supertype as the standard potential, hence, they share most of their properties and methods.

potential_dft.show(cbar=True);

Comparing DFT to IAM

To show the difference that including charge trasnfer can have, we compare the DFT potential to an equivalent potential using the IAM. We make sure to set projection="finite" to ensure that the projection integrals are done identically.

atoms = gpaw.atoms

potential_iam = (
    abtem.Potential(atoms, gpts=potential_dft.gpts, projection="finite")
    .build()
    .compute()
)

Below compare the projection of the potentials. It is difficult to discern any visual differences between the IAM and DFT model from their resepective heatmaps, hence, we also calculate and show their difference and relative difference. In hBN the negative charge is transferred from Boron to Nitrogen, thus better screening the Nitrogen cores and lowering the potential.

To make the comparison easier, we subtract the minimum from both potentials, this does does change the physics.

projected_potential_dft = potential_dft.project()
projected_potential_iam = potential_iam.project()

projected_potential_dft -= projected_potential_dft.min()
projected_potential_iam -= projected_potential_iam.min()

difference = projected_potential_dft - projected_potential_iam

relative_difference = projected_potential_dft.relative_difference(
    projected_potential_iam, min_relative_tol=0.001
)

fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize=(13, 4), sharey=True)

projected_potential_iam.show(ax=ax1, cbar=True, title="iam potential", vmax=400)

projected_potential_dft.show(ax=ax2, cbar=True, title="dft potential", vmax=400)

difference.show(
    ax=ax3,
    cbar=True,
    title="difference (dft - iam)",
    vmin=-12,
    cmap="seismic",
    vmax=12,
)

relative_difference.show(
    ax=ax4,
    cbar=True,
    title="difference (dft - iam)",
    vmin=-60,
    cmap="seismic",
    vmax=60,
)

plt.tight_layout()

Comparing lineprofiles as below is often a preferable way of showing differences between two images. Below we show the lineprofiles of \(x=0\).

iam_line = projected_potential_iam.interpolate_line(
    start=(0, 0), end=(0, projected_potential_iam.extent[1])
)

dft_line = projected_potential_dft.interpolate_line(
    start=(0, 0), end=(0, projected_potential_dft.extent[1])
)

abtem.stack([iam_line, dft_line], ("IAM", "DFT")).show();

Electron diffraction of DFT potentials

Next, we calculate the influence of charge transfer on the electron diffraction patterns in hBN, approximately reproducing a published result [SML+19].

We create a PlaneWave at an energy of \(100 \ \mathrm{keV}\), then run the multislice algorithm for the IAM and DFT potential, finally caluclating the diffraction patterns.

planewave = abtem.PlaneWave(energy=100e3)

diffraction_iam = (
    planewave.multislice(potential_iam)
    .diffraction_patterns(block_direct=True, max_angle=50)
    .compute()
)
diffraction_dft = (
    planewave.multislice(potential_dft)
    .diffraction_patterns(block_direct=True, max_angle=50)
    .compute()
)

stacked = abtem.stack([diffraction_iam, diffraction_dft], ("iam", "dft"))

stacked.show(explode=True, common_color_scale=True, cbar=True, figsize=(12, 6));

We show the resulting diffraction patterns below. The most obvious difference is that order of intensity of the \((200)\) and \((110)\) families of diffraction spots is flipped in the two models.

diffraction_spots = stacked.index_diffraction_spots(atoms)

visualization = diffraction_spots.show(
    explode=True,
    common_color_scale=False,
    cbar=True,
    figsize=(13, 6),
    scale=0.02,
    annotation_kwargs={"threshold": 0.000001},
)

Below we compare the relative intensities of the diffraction spots in a table, the intensities are normalized to the intensity of the \((110)\) spots. We also include experimental values cite{gpawpotential}.

We see the largest difference between the models for the \((200)\) spots. Comparing to the experimental result, we see that IAM predicts the wrong intensity ordering of the \((110)\) and \((200)\) spots. The DFT model gets the order of the spots right, however, the relative intensity is significantly overestimated. The main cause of the discrepancy is the noninclusion of phonon scattering, as we shall see in the following section.

diffraction_spots.to_dataframe()

hkl	0 0 0	0 1 0	0 2 0	0 3 0	0 -3 0	0 -2 0	0 -1 0	1 0 0	1 1 0	1 2 0	...	-2 1 0	-2 2 0	-2 3 0	-2 -1 0	-1 0 0	-1 1 0	-1 2 0	-1 3 0	-1 -2 0	-1 -1 0
iam	0.0	0.000106	0.000016	0.000014	0.000014	0.000016	0.000100	0.000100	0.000094	0.000006	...	0.000094	0.000016	0.000006	0.000006	0.000106	0.000100	0.000094	0.000006	0.000006	0.000094
dft	0.0	0.000088	0.000017	0.000014	0.000014	0.000018	0.000083	0.000083	0.000096	0.000006	...	0.000096	0.000018	0.000006	0.000006	0.000088	0.000083	0.000096	0.000006	0.000006	0.000096

2 rows × 37 columns

Frozen phonons and DFT

Currently Under Construction

Check back for more…