Semi Empirical

Introduction

SemiEmpirical (SE) calculator can model the electronic properties of molecules, crystals and devices using both self-consistent and non-self-consistent tight-binding models. In this chapter, the implemented tight-binding models based on the Slater–Koster model and the extended Hückel model are presented.

The Slater–Koster tight-binding model follows closely the DFTB formalism described in [1], and it is recommended that this paper, and [2], are cited in publications using the SemiEmpiricalCalculator and DeviceSemiEmpiricalCalculator with the SlaterKosterHamiltonianParametrization in QuantumATK.

The extended Hückel model in SE is described in [2], and it is recommended that this paper is cited in publications using the SemiEmpiricalCalculator and DeviceSemiEmpiricalCalculator with the HuckelHamiltonianParametrization in QuantumATK.

In SE, the non-self-consistent part of the tight-binding Hamiltonian is parametrized using a two-center approximation, i.e. the matrix elements only depend on the distance between two atoms and is independent of the position of the other atoms. In the extended Hückel model, the matrix elements are described in terms of overlaps between Slater orbitals on each site. In this way, the matrix elements can be defined by very few parameters. In the Slater–Koster model, the distance-dependence of the matrix elements is given as a numerical function; this gives higher flexibility but also makes the fitting procedure more difficult.

The self-consistent part of the calculation is identical for both SE models. The density matrix is calculated from the Hamiltonian using non-equilibrium Green’s functions for device systems, while for molecules and crystals it is calculated by diagonalization. The density matrix defines the real-space electron density and consequently the Hartree potential can be obtained by solving the Poisson equation. The following describes the details of the mathematical formalism behind the implementation.

For a list of the built-in parameter sets in SE, see Built-in parameter sets in ATK-SE. In addition to these, it is possible to directly use parameter sets downloaded from the DFTB website or define your own Slater–Koster table or extended Hückel basis set in the QuantumATK input scripts.

Background information

Non-self-consistent Hamiltonian

The Hamiltonian is expanded in a basis of local atomic orbitals (an LCAO expansion),

\[\phi_{nlm}(\mathbf{r}) = R_{nl}(r) Y_{lm}(\hat{\mathbf{r}}),\]

where \(Y_{lm}\) is a spherical harmonic and \(R_{nl}\) is a radial function. Typically, the atomic orbitals used in the LCAO expansion have a close resemblance to the atomic eigenfunctions.

Onsite terms

With this form of the basis set, the onsite elements are given by

\[S_{ij}^\mathrm{onsite} = \delta_{ij},\]
\[H_{ij}^\mathrm{onsite} = E_i \delta_{ij},\]

where \(E_i\) is an adjustable parameter that is often close to the atomic eigenenergy.

Offsite terms in the extended Hückel model

The central object in the extended Hückel model is the overlap matrix:

\[S_{ij} = \int_V \phi_i(\mathbf{r} -\mathbf{R}_i) \phi_j(\mathbf{r}- \mathbf{R}_j) \; \mathrm{d}\mathbf{r}.\]

To calculate this integral, the form of the basis functions must be specified. In the extended Hückel model, the basis functions are parametrized by Slater orbitals,

\[R_{nl}(r) = \frac{r^{n-1}}{\sqrt{(2n)!}} \left[C_1 (2 \eta_1)^{n+\frac{1}{2}} e^{-\eta_1 r}+C_2 (2 \eta_2)^{n+\frac{1}{2}} e^{-\eta_2 r} \right].\]

The LCAO basis is described by the adjustable parameters \(\eta_1\), \(\eta_2\), \(C_1\), and \(C_2\). These parameters must be defined for each angular shell of valence orbitals for each element.

The overlap matrix defines the Hamiltonian,

\[H_{ij}= \frac{1}{4} (\beta_i+\beta_j) (E_i+E_j) S_{ij},\]

where \(E_i\) is the onsite orbital energy and \(\beta_i\) is a Hückel fitting parameter (often chosen to be 1.75).

Weighting schemes

The are two variants of the weighting schemes for the orbital energies of the offsite Hamiltonian. The scheme used above,

\[\frac{1}{2} \beta (E_i + E_j) S_{ij},\]

where \(\beta = \frac{1}{2} (\beta_i+\beta_j)\) is used by Wolfsberg [3], while Hoffmann ([4] and [5]) uses

\[\frac{1}{2} (\beta +\alpha^2+(1-\beta) \alpha^4) (E_i + E_j),\]

where \(\alpha = (E_i - E_j)/(E_i + E_j)\).

Both variants are available in Semi Empirical through the parameters WolfsbergWeighting and HoffmannWeighting, which can be given to the HuckelHamiltonianParametrization class.

Offsite Hamiltonian in the Slater–Koster model

The overlap matrix is given by pairwise integrals between the different basis functions. These integrals can be pre-calculated for all relevant distances and different orbital combinations, and stored in so-called Slater–Koster tables. The Slater–Koster table stores the distance-dependent parameters \(s(d, Z_1, Z_2, l_1, l_2, m)\), where \(d\) is the distance, \(Z_1\), \(Z_2\) the element types, \(l_1\), \(l_2\) the angular momentum of the two orbitals, and the index \(m \le \min(l_1, l_2)\).

From the Slater–Koster tables, the overlap matrix elements are given by

\[S_{ij} = \sum_{m \le \min(l_i, l_j)} \alpha_{l_i, m_i, l_j, m_j, m}(\hat{R}_{ij}) s(d_{ij}, Z_i, Z_j, l_i, l_j, m),\]

where \(\alpha\) are the Slater–Koster expansion coefficients.

In the Slater–Koster, model it is assumed that also the Hamiltonian has a pairwise form and a Slater–Koster table is generated for the Hamiltonian matrix elements. This table may be generated by calculating Hamiltonian matrix elements for a set of dimer distances or by simply fitting matrix elements to the band structure for different lattice constants.

In Semi Empirical, the Slater–Koster table is constructed either by providing the path to a directory containing compatible Slater–Koster files (see DFTBDirectory and HotbitDirectory), or directly using the SlaterKosterTable class.

Note that the extended Hückel model is a Slater–Koster model too, but with a special fitting procedure for the Hamiltonian matrix elements.

Self-consistent Hamiltonian

In the self-consistent semi-empirical models in QuantumATK, the electron density is computed using the tight-binding model as described above. The density gives rise to a Hartree potential \(V_H\). The calculation of the Hartree potential is described in detail in the section on Poisson solvers.

The Hartree potential is included through an additional term in the Hamiltonian:

\[H_{ij}^\mathrm{SCF} = \frac{1}{2}(V_H(\mathbf{R}_i)+V_H(\mathbf{R}_j)) S_{ij}.\]

Electron density

The electron density is given by the occupied eigenfunctions,

\[n(\mathbf{r}) = \sum_{\alpha} f_\alpha |\psi_\alpha(\mathbf{r})|^2,\]

where \(f_\alpha\) is the occupation of the level denoted by \(\alpha\). For finite temperature calculation the occupations are determined by the Fermi-Dirac distribution \(f_\alpha = \frac{1}{1 + e^{(\epsilon_\alpha - \epsilon_F)/kT}}\) with \(\epsilon_\alpha\) being the energy of the eigenstate \(\psi_\alpha\), \(\epsilon_F\) the Fermi level and \(T\) the electron temperature. However, other smooth distributions may be introduced in order to help speed convergence (see Occupation Methods).

The eigenstates in the Slater orbital basis can be written as

\[\psi_\alpha = \sum_i c_{\alpha i} \phi_i,\]

where the total number of electrons, \(N=\int_V n(\mathbf{r}) \, \mathrm{d}\mathbf{r}\), is given by

\[N = \sum_{i j} D_{ij} S_{ij},\]

where \(D_{ij} = \sum_{\alpha} f_\alpha c_{\alpha i}^* c_{\alpha j}\) is the density matrix.

An approximate atom-based electron density

In practice, a simple approximation is used for the electron density. To this end, we introduce the Mulliken population,

\[m_l = \sum_{i \in l} \sum_j D_{ij} S_{ij},\]

for shell \(l\) of atom number \(\mu\), and write the total number of electrons as a sum of atomic contributions, \(N=\sum_\mu \sum_{l \in \mu} m_l\). The radial dependence of each atomic-like density is represented by a Gaussian function, and the total induced charge in the system is approximated by

\[\delta n(\mathbf{r})=\sum_\mu \sum_{l \in \mu} \delta m_l \left({ \frac{\alpha_l}{\pi}}\right)^{\frac{3}{2}} e^{-\alpha_l |\mathbf{r}-\mathbf{R}_\mu|^2},\]

where \(\delta m_l = m_l - Z_\mu\) is the total charge for shell \(l\) of atom \(\mu\), i.e. the sum of the valence electron charge \(m_l\) and the ionic charge \(-Z_\mu\).

To see the significance of the width \(\alpha_l\) of the Gaussian orbital, consider the electrostatic potential from a single Gaussian density at position \(\mathbf{R}_\mu\):

\[V_{H}(\mathbf{r}) = (m_l-Z_\mu) \frac{\text{erf}(\sqrt{\alpha_l} | \mathbf{r}-\mathbf{R}_\mu|) }{ | \mathbf{r}-\mathbf{R}_\mu| }.\]

The onsite value of the Hartree potential is \(V_{H}(\mathbf{R}_\mu)=(m_l-Z_\mu) U_l\), where \(U_l= 2 \sqrt{ \frac{\alpha_l}{\pi}}\) is the onsite Hartree shift. In Semi Empirical, it is the specified value of \(U_l\) that is used to determine the width \(\alpha_l\) of the Gaussian using the above relation.

Onsite Hartree shift parameters

The shell-dependent onsite Hartree shift \(U_l\) can be obtained from an atomic calculation. \(U_l\) is related to the linear shift of the eigenenergy \(\varepsilon_l\) of shell \(l\) as function of the shell occupation \(q_l\):

\[U_l = \frac{d \varepsilon_l}{d q_l}.\]

Thus, \(U_l\) can be obtained by performing atomic calculations with different values of \(q_l\).

In QuantumATK, it is recommended to use the same onsite Hartree shift parameter for the s- and p-shells of each atom.

ATK provides a database for \(U_l\) calculated using DFT and the PBE functional (see Generalized-gradient approximation (GGA)). Access to the data is through the function ATK_U. Due to backwards compatibility, the HoffmanHückelParameters and MullerHückelParameters do not use the ATK_U database, but use special values of the electrostatic parameter \(U\), see Element data.

Spin polarization

The inclusion of spin in the tight-binding Hamiltonian follows the scheme in [6]. The following spin dependent term is added to the Hamiltonian:

\[H_{ij}^\sigma= \pm \frac{1}{2} S_{ij} \left(dE_{l_i} + dE_{l_j} \right),\]

where the sign in the equation depends on the spin.

The spin splitting \(dE_{l}\) of shell \(l\) is calculated from the spin-dependent Mulliken populations \(\mu_l\) of each shell at the local site:

\[dE_{l} = \sum_{l' \in \mu_l} W_{l l'} \, (m_{l'\uparrow}- m_{l'\downarrow}).\]

Onsite spin-split parameters

The shell-dependent spin splitting strength \(W_{ll'}\) can be obtained from a spin-polarized atomic calculation [6]:

\[W_{l l'} = \frac{1}{2}\left(\frac{d \varepsilon_{l\uparrow}}{d m_{l' \uparrow}} -\frac{d \varepsilon_{l\uparrow}}{d m_{l'\downarrow}} \right).\]

Since \(W_{l l'}\) enters symmetrically in the Hamiltonian, it is convenient to symmetrize it:

\[\bar{W}_{l l'} = \frac{1}{2} (W_{l l'} + W_{l' l}).\]

ATK provides a database for \(\bar{W}_{l l'}\). Access to the data is through the function ATK_W.

Tight-binding total energy

The calculation of the total energy follows [7] and [1]. The total energy has five terms:

\[E=E_{\rm H^0}+E_{\rm\delta H}+E_{\rm ext}+E_{\rm spin} +E_{\rm pp}.\]
  • \(E_{\rm H^0}\) is the one-electron energy of the non-self-consistent Hamiltonian, given by

    \[E_{\rm H^0} = \sum_{ij} D_{ij} H^0_{ij}.\]
  • \(E_{\rm\delta H}\) is the electrostatic difference energy,

    \[E_{\rm\delta H} = \frac{1}{2}\int V^H_0(\mathbf{r}) \delta n(\mathbf{r})d\mathbf{r}.\]
  • \(E_{\rm ext}\) is the electrostatic interaction between the electrons and an external field,

    \[E_{\rm ext} = \int V^{\rm ext}(\mathbf{r}) \delta n(\mathbf{r})d\mathbf{r}.\]
  • \(E_{\rm spin}\) is the spin polarization energy [6],

    \[E_{\rm spin} = -\frac{1}{2} \sum_\mu \sum_{l \in \mu} \sum_{l' \in \mu} W(Z_{\mu}, l,l') m_l m_{l'}.\]
  • \(E_{\rm pp}\) is the repulsive energy from a pair-potential between each atom pair, \(V^{\rm pp}(Z_{\mu}, Z_{\mu'},R_{\mu, \mu'})\),

    \[E_{\rm pp} = \sum_{\mu < \mu'} V^{\rm pp}(Z_{\mu}, Z_{\mu'},R_{\mu, \mu'}).\]

    It is optional to add this term to the tight-binding model; it does not affect the electronic structure. The tight-binding model will, however, not give sensible geometries without a repulsive pair-potential.

Important

In the current version of QuantumATK, only the DFTB and Hotbit parameter sets contain a repulsive pair potential term, and only these methods can be used for geometry optimizations.

Parameters

Parameters for the Slater–Koster method

The Slater–Koster parameters can be provided either through the SlaterKosterTable class or through various 3rd-party formats. The supported 3rd-party formats are the DFTB Slater–Koster files from the DFTB consortium, and the Hotbit Slater–Koster files from the Hotbit consortium.

To use parameters in the DFTB or Hotbit format, put the parameter files in a single directory, and setup the basis set and pair potentials using the functions DFTBDirectory or HotbitDirectory.

Tables showing which elements are covered by the parameters are given here: Slater–Koster basis sets.

Shipped DFTB and Hotbit parameters

The current version of SE is shipped with a number of DFTB style parameters from the CP2K and Hotbit consortia. We recommended that these resources are cited if the parameters are used in publications. It is most easy to setup these basis sets using the QuantumATK interface.

Shipped Slater–Koster Table parameters

A number of orthogonal tight-binding parameters are provided. The parameters are from Vogl et al. [8] and Jancu et al. [9], and it is recommended that these papers are cited if the parameters are used for publications. These basis sets can easily be set up using QuantumATK.

Parameters for the extended Hückel method

The parameters \(\eta_1\), \(\eta_2\), \(C_1\), \(C_2\), and \(E\) must be defined for each valence orbital, while \(\beta\) and \(U\) only depend on the element type. Different parameter sets are provided with Semi Empirical, but it is also possible to provide user-defined parameters in the input file using the HuckelBasisParameters class.

The tables below provide a mapping between the symbols in the equations and the corresponding keywords.

Table 12 HuckelBasisParameters.

Symbol

HuckelBasisParameters

\(E_i\)

ionization_potential

\(\beta\)

wolfsberg_helmholtz_constant

\(U\)

onsite_hartree_shift

\(W\)

onsite_spin_split

\(E^{\mathrm{VAC}}\)

vacuum_level


Table 13 SlaterOrbital parameters.

Symbol

SlaterOrbital parameters

\(n\)

principal_quantum_number

\(l\)

angular_momentum

\(\eta\)

slater_coefficients

\(C\)

weights

The current version of QuantumATK comes with built-in Hoffmann and Müller parameter sets, which are appropriate for organic molecules. For crystalline structures, both metals and organic materials like graphene, parameters from J. Cerda are provided. When using these parameters, [10] should be referenced. Tables with the available parameter sets can be found here: Extended Hückel basis sets.

To combine parameters from different sources, it is important to make sure they use the same energy zero level, in order to obtain correct charge transfers. This can be obtained by ensuring that the crystals have the correct work function and molecules the correct ionisation energies. For this purpose, an additional parameter \(E^{\mathrm{VAC}}\) is introduced, which shifts the energy of the vacuum level: If a calculation with \(E^{\mathrm{VAC}}=0\) eV has a work function of 6.5 eV, then by setting \(E^{\mathrm{VAC}}=-1.5\) eV all bands shift rigidly upwards by 1.5 eV, and the work function becomes 5.0 eV.

Note

The Hückel parameters have been fitted for non-self-consistent calculations. To use the parameters in self-consistent calculations, the self-consistent onsite shifts must be compensated by a reverse shift of the vacuum_levels.