The NRL Tight-Binding Codes

---------------------------------------

The NRL Tight-binding method provides an efficient method for calculating properties of materials. The advantage of the NRL-TB method over classical potential simulations is that it explicitly incorporates the real electronic structure and bonding of the material, obtained by an interpolation from a database of first-principles results.

There are now thirteen codes used to perform these calculations:

  1. Fitting code spd for determining tight-binding parameters from the first-principles database. This code has been shown to produce accurate tight-binding parameters for elemental materials and is extended to binary materials.
  2. static, a total energy and electronic structure evaluation code which can be used for metals, semi-conductors and insulators, but does not do molecular dynamics.
  3. TBMD, a Tight Binding Molecular Dynamics code for metals. Developed by Florian Kirchhoff.
  4. SCTB, a charge self-consistent tight-binding total energy evaluation code.
  5. Fitting code extension to treat f-states.
  6. Static code extension to treat f-states.
  7. Static-Spin-Orbit f-states
  8. Spin-Polarized-Fit spd
  9. Spin-Polarized-Static spd
  10. Coherent Potential Approximation
  11. Cluster Energies
  12. Spin-Orbit (p and d levels)
  13. Phonon Calculations

---------------------------------------

Theory of the NRL Tight-Binding Method

The NRL-Tight-Binding method (TB) is a Slater-Koster-like TB method [[1]–[3]], which has been successfully applied to many single element, binary and ternary systems [[4], [5]]. The method performs very well in metals, insulators, and semiconductors [[6], [7]], and can be extended to clusters and molecules. It is built on a fitting process to Density Functional Theory results. It works in both an orthogonal and a nonorthogonal basis and computes both the energy bands and total energies. It has transferability far beyond the fitted DFT database. The NRL-TB formalism with the basic equations is outlined below. First, the total energy is determined using the fact that in DFT the total energy $E[n(r)]$ is given by the expression:


\begin{displaymath}E[n(r)] = \sum_{i} \epsilon_{i} + F[n(r)] \end{displaymath}

where $n(r)$ is the electron density and $F[n(r)]$ contains the remaining parts of the DFT total energy minus the sum of the one-electron eigenvalues. In other TB approaches $F[n(r)]$ is a sum of pair potentials. However, based on the fact that in DFT one can uniformly shift the energy bands by a constant, in this method a shift $V_{0} = F[n(r)]/N$ is applied to the first-principles eigenvalues $\epsilon_{i}$, where $N$ is the number of valence electrons in the system. This shifts the one-electron eigenvalues to new values,


\begin{displaymath}\epsilon_{i}^{\prime} = \epsilon_{i} + V_{0} \end{displaymath}

and hence the total energy becomes the sum of the shifted DFT eigenvalues over the Brillouin zone, i.e.


\begin{displaymath}E[n(r)] = \sum_{i} \epsilon_{i}^{\prime} \end{displaymath}

This is a two-center TB scheme, where the on-site terms have a polynomial form as a function of the atomic density. For a single element, the density of atom i is defined as


\begin{displaymath}\rho_{i} = \sum_{i} \exp( -\lambda^2 R_{ij}) F_{C}(R_{ij}) \end{displaymath}

where the sum is over all the neighboring atoms j within a range of cutoff distance $R_{c}$ of atom $i$, is a fitting parameter, and $F_{C}(R_{ij})$ is a smooth cutoff function. The angular-momentum-dependent on-site terms are defined by


\begin{displaymath}h_{\ell}(R_{i}) = \alpha_{i\ell} + \beta_{i\ell}\rho^{2/3} + \gamma_{i\ell}\rho^{4/3} + \delta_{i\ell}\rho^{2} \end{displaymath}

where $\ell$ represents the $s$, $p$, and $d$ orbitals, and $\alpha\ell$, $\beta\ell$, $\gamma\ell$ and $\delta\ell$ are fitting coefficients. The two-center $s-p-d$ Slater-Koster (SK) Hamiltonian and overlap integrals consist of ten independent SK parameters, which are assumed to all have polynomial times exponential forms in terms of the neighbor distance R, given by


\begin{displaymath}P_{\gamma}(R) = (e_{\gamma} + f_{\gamma}R + g_{\gamma}R^{2}) \exp([-q_{\gamma}^{2}]F_{C}(R)) \end{displaymath}

where $\gamma$ indicates the type of ineractions, including $ss\sigma$, $pp\sigma$, $sp\sigma$, $dd\sigma$, $sd\sigma$, $pd\sigma$, $pp\pi$, $dd\pi$, $pd\pi$, and $dd\delta$. $e_{\gamma}$, $f_{\gamma}$, $g_{\gamma}$, and $q_{\gamma}$ are our fitting coefficients. The SK overlap functions, in a non-orthogonal calculation, are defined to have the same form as the hopping parameters. Overall, there are in total 93 fitting coefficients for a single element in the on-site, hopping and overlap terms in the TB Hamiltonian with s, p, and d orbitals.

For semi-conductors, in the overlap matrix we use a polynomial of the form


\begin{displaymath}S_{\gamma}(R) = ( \delta_{\gamma} + e_{\gamma} + f_{\gamma}R + g_{\gamma}R^{2}) \exp([-q_{\gamma}^{2}]F_{C}(R)) \end{displaymath}

where $\delta_{\gamma}$ is the Kronecker delta.

In order to determine the above coefficients a least-squares procedure is used to fit to DFT total energies and energy bands as a function of volume for high symmetry structures as bcc, fcc, and simple cubic crystal structures with varying atomic volumes. In some cases, to improve the fit one uses selected volumes for the hcp structure, and for semiconductors fit the diamond lattice as well. The total energy is usually weighed at around 200-300 times over a single band energy. In general, the fitting RMS error is less than 10 mRy and 0.2 mRy for the energy bands and total energy, respectively.

Having determined the above coefficients, the method is used to predict total energies accessible to standard DFT which were not fitted, such as elastic constants, phonon spectra and surface energies. Furthermore, large scale simulations can be performed which are not practical via DFT, such as static calculations for the energetics of systems containing up to 20,000 atoms, or calculations for very large number of k-points needed for mapping Fermi surfaces and evaluating susceptibilities.

The method also has capabilities to perform molecular dynamics (MD) simulations accommodating as many as 10,000 atoms, and 10,000 MD steps, which is an impossible task for standard DFT codes. Such calculations should yield mean square displacements, thermal expansion, and vacancy formation energies [[2]–[4]].

References

1
J. C. Slater and G. F. Koster, Phys. Rev. 94, 1498 (1954).

2
R. E. Cohen, M. J. Mehl, and D. A. Papaconstantopoulos, Phys. Rev. B 50, 14694 (1994).

3
M. J. Mehl, and D. A. Papaconstantopoulos, Phys. Rev. B 54, 4519 (1996).

4
D. A. Papaconstantopoulos, and M. J. Mehl, J. Phys.-Condes. Matter 15, R413 (2003).

5
D. A. Papaconstantopoulos, M. J. Mehl, and M. D. Johannes, Phys. Rev. B 82, 054503 (2010).

6
N. Bernstein, M. J. Mehl, D. A. Papaconstantopoulos, N. I. Papanicolaou, M. Z. Bazant, and E. Kaxiras, Phys. Rev. B 62, 4477 (2000).

7
N. Bernstein, M. J. Mehl, and D. A. Papaconstantopoulos, Phys. Rev. B 66, 075212 (2002).

8
S. Silayi, D. A. Papaconstantopoulos, and M.J. Mehl Comput. Mater. Sci. 146, 278–286 (2018)

9
S. Silayi, P-H. Chang, and D. A. Papaconstantopoulos, Comput. Mater. Sci. 173, 109454 (2020)