Removing all periodic boundary conditions: Efficient non-equilibrium Green function calculations

TL;DR

This paper introduces a computational method that calculates electronic structure and transport without relying on periodic boundary conditions, reducing artifacts and computational costs in non-equilibrium Green function simulations.

Contribution

The authors develop and implement a new approach that eliminates the need for periodic super-cells in Green function calculations, improving accuracy and efficiency.

Findings

Reduces CPU-hours compared to super-cell methods

Successfully applied to graphene defect and junction examples

Eliminates spurious periodic images and interference effects

Abstract

We describe a method and its implementation for calculating electronic structure and electron transport without approximating the structure using periodic super-cells. This effectively removes spurious periodic images and interference effects. Our method is based on already established methods readily available in the non-equilibrium Green function formalism and allows for non-equilibrium transport. We present examples of a N defect in graphene, finite voltage bias transport in a point-contact to graphene, and a graphene-nanoribbon junction. This method is less costly, in terms of CPU-hours, than the super-cell approximation.

Tables

Table 1. Table 1: Maximum disk space requirements for the 3 test examples. TS: TranSiesta ; TBT: TBtrans . For 𝐤 𝐤 \mathbf{k} resolved TBT self-energy calculations one need not store the self-energies for the since those can be calculated when needed. For 300 300 300 𝐤 𝐤 \mathbf{k} points they would use ∼ 500 ​ GB similar-to absent 500 GB \sim 500\,\mathrm{GB} .

Graphene-Nitrogen	TS	TBT
$n$, $n_{𝐤}$	$n_E=51$	$n_E=600$
${𝚺}^{ℛ}$: $648$, $6$	$1.99\mathrm{GB}$	—
${𝚺}^{ℛ}$: $2052$, $1$	$3.33\mathrm{GB}$	$37.8\mathrm{GB}$
Graphene-STM	TS	TBT
$n$, $n_{𝐤}$	$n_E=178$	$n_E=600$
${𝚺}^{ℛ}$: $720$, $6$	$8.34\mathrm{GB}$	—
${𝚺}^{ℛ}$: $2160$, $1$	$12.5\mathrm{GB}$	$41.9\mathrm{GB}$
Graphene-zGNR	TS	TBT
$n$, $n_{𝐤}$	$n_E=178$	$n_E=600$
${𝚺}^{ℛ}$: $648$, $6$	$6.76\mathrm{GB}$	—
${𝚺}^{ℛ}$: $1060$, $1$	$3.01\mathrm{GB}$	$10.1\mathrm{GB}$

Table 2. Table 2: Timings of the various steps in the presented calculations, the timings are seconds per core in a 20 core setup. TS: TranSiesta ; TBT: TBtrans ; 𝚺 ℛ superscript 𝚺 ℛ \boldsymbol{\Sigma}^{\mathcal{R}} : calculating the real-space self-energies for both TS and TBT. All calculations are done on the same machine.

Graphene-Nitrogen	Timing [s]	Total [s]
Siesta  $+$ PDOS	$19187$	$19187$
TS $+$ TBT	$1885+5484$	$7369$
TS $+$ TBT $+$ ${𝚺}^{ℛ}$	$1364+324+1410$	$3099$
Graphene-STM	Timing [s]	Total [s]
TS $+$ TBT	$5372+31043$	$36415$
TS $+$ TBT $+$ ${𝚺}^{ℛ}$	$5465+248+3608$	$9322$
Graphene-zGNR	Timing [s]	Total [s]
TS $+$ TBT	$10180+8206$	$18386$
TS $+$ TBT $+$ ${𝚺}^{ℛ}$	$2836+157+7409$	$10401$

Equations

$${\mathbf{\Sigma }}_{00}^{\mathcal{R}}={\mathbf{S}}^{\mathcal{R}}z-{\mathbf{H}}^{\mathcal{R}}-[{\mathbf{G}}_{00}^{\mathcal{R}}{]}^{-1}.$$

$${\mathbf{H}}^{\mathcal{R}}=\left[\begin{array}{cccc}{\mathbf{H}}_{00}& {\mathbf{H}}_{01}& 0& \dots \\ {\mathbf{H}}_{10}& {\mathbf{H}}_{11}& {\mathbf{H}}_{12}& \\ 0& {\mathbf{H}}_{21}& {\mathbf{H}}_{22}& \\ ⋮& & & \ddots \end{array}\right].$$

$${\mathbf{G}}_{00}^{\mathcal{R}}(z)$$

$$=[{\mathbf{S}}^{\mathcal{R}}z-{\mathbf{H}}^{\mathcal{R}}{]}_{00}^{-1}=[{\mathbf{S}}_{00}z-{\mathbf{H}}_{00}-{\mathbf{\Sigma }}_{00}^{\mathcal{R}}{]}^{-1}$$

$${\mathbf{\Sigma }}_0^L$$

$${\mathbf{L}}_0$$

$${\mathbf{R}}_0$$

$${\mathbf{\tau }}^L$$

$${\mathbf{\tau }}^R$$

$${\mathbf{\Sigma }}_i^L$$

$${\mathbf{\Sigma }}_i^R$$

$${\mathbf{L}}_i$$

$${\mathbf{R}}_i$$

$$\mathbf{H}=\left[\begin{array}{cccc}{\mathbf{H}}_{00}& {\mathbf{H}}_{01}& 0& \\ {\mathbf{H}}_{10}& {\mathbf{H}}_{00}& {\mathbf{H}}_{01}& \\ 0& {\mathbf{H}}_{10}& \ddots & \end{array}\right].$$

$$\tilde{\mathbf{Y}}$$

$$\tilde{\mathbf{X}}$$

$${\mathbf{G}}_{nn}$$

$${\mathbf{G}}_{mn}$$

$${\mathbf{G}}_K^N=\frac{1}{N}\sum _{\begin{array}{c}j=0\\ k_j=\frac{K+j}{N}\end{array}}^{N-1}\left[\begin{array}{ccc}1& \cdots & e^{i(1-N)k_j}\\ e^{i{k}_j}& \cdots & e^{i(2-N)k_j}\\ ⋮& \ddots & ⋮\\ e^{i(N-1)k_j}& \cdots & 1\end{array}\right]\otimes {\mathbf{G}}_{k_j}^1$$

$${\mathbf{\Sigma }}_{00}^{\mathcal{R}}$$

$$M_{\mathrm{GB}}(n,n_{\mathbf{k}},n_E)=16\cdot n_{\mathbf{k}}{n}^2(2+n_E)/1024^3,$$

Full text

I Introduction

Widely used, efficient computational methods have been developed for calculations of electronic structure of systems presenting perfect periodic repetition of a unit cell along one, two, or three dimensions surrounded by vacuum. These are typically based on Kohn-Sham Density Functional Theory (DFT)Martin (2004); Kohanoff (2006). The infinite system is replaced by finite unit-cell with periodic boundary conditions (PBC) using Bloch’s theorem and a discrete sampling of Bloch phases/$\mathbf{k}$-points. Due to the efficient implementations this method is also applied systems which lack periodicity. For example surfaces are modelled by a slab, isolated defects by periodically repeated defects surrounded by “large” regions of bulk, and isolated adsorbates on surfaces by a mix. This results in compromises due to computational feasibility with respect to slab-size and inter-defect distances, which may lead to unwanted effects related to interferences or standing-wave patterns not present in the ideal, large system. Beyond PBC methods have been around for a long time. These include matching of the wavefunctions in different regions, e.g. surface and bulkAppelbaum and Hamann (1972), and Green function or embedding methodsWilliams et al. (1982); MacLaren et al. (1989); Feibelman (1992); Ishida (1997); Inglesfield (2001); Li et al. (2018) have e.g. been used to treat the isolated defect/adsorbate on a surface or electronic transport between two electrodesWortmann et al. (2002). These methods are based on a screening assumption where the potential has converged to its bulk value outside the computational “active” region.

In particular, for transport calculations the treatment of systems as “open” with semi-infinite electrodes along the transport directions is essential. A number of computational implementations have been developed for more than a decade for this problem based on the non-equilibrium Green function (NEGF) methodTaylor et al. (2001); Brandbyge et al. (2002); Thygesen and Jacobsen (2005); Papior et al. (2017); Ozaki et al. (2010); Smidstrup et al. (2019). These typically represent electrodes by a unit-cell repeated as periodic layers along the semi-infinite electrode/transport direction, and use PBC and corresponding $\mathbf{k}$-points in the directions transverse to this. A self-energy is then used to treat the semi-infinite direction in a numerically exact way based on a very efficient methodSancho et al. (1984) which recursively removes the infinite number of degrees of freedom/states in the semi-infinite direction. This approach may also be used in “single-electrode mode” treating the surface of semi-infinite bulk with a computational load comparable to slab calculation of e.g. chemical reactions at the surface Sánchez-Portal (2007); Papior et al. (2017); Smidstrup et al. (2017). Indeed this avoids the periodic images and finite size effects of the slabs in the surface-normal direction, but leaves the periodicity in the surface direction.

In this paper we present a simple, efficient and precise method based on Green function theory which can be used for isolated defects as well as extended NEGF calculations using multiple probes/electrodes. Our method solves this problem by calculating the real-space self-energy which can be outlined as (details explained in the Method section)

[TABLE]

Equation (1) is computationally demanding since real-space quantities requires a dense integration grid in reciprocal space. We emphasize that our method focuses on the efficient algorithmic implementation which has prohibited the community to extend its broader use. Secondly, our method allows non-equilibrium calculations by the regular assumption of “equilibrium” electrodesTaylor et al. (2001); Brandbyge et al. (2002); Thygesen and Jacobsen (2005); Papior et al. (2017); Ozaki et al. (2010); Smidstrup et al. (2019). The paper is organized as follows. First we describe the theoretical and computational details of our method. Then we show DFT$+$NEGF results using the real-space self-energies for three illustrative cases: i) electronic structure of a nitrogen defect in a large graphene lattice, ii) non-equilibrium transport in a gold STM tip in contact with a graphene flake, and iii) a graphene/graphene nanoribbon junction.

II Method

We remark that Eq. (1) is a well known equation in the transport community and that our contribution here lies in the implementation. In the following we will describe the method for a pristine bulk system which is the basis for defected systems.

The starting point of the efficient real-space self-energy method is any system with PBC in 2 or 3 directions where one wishes to replace a predefined direction with a semi infinite description, see Fig. 1a. This may be efficiently described using two semi-infinite directions and one PBC direction, see Fig. 1b, lastly our presented method replaces any number of PBC, and/or semi-infinite directions, with a single self-energy, see Fig. 1c.

The Hamiltonian of an infinitely large system may be written in a block-tri-diagonal “shell” fashion:

[TABLE]

Here index [math] in $\mathbf{H}_{00}$ is referred as the primary unit cell with only nearest neighbour couplings, $\mathbf{H}_{01}$ is the coupling between the primary unit-cell and the first set of neighbour cells (2D: $8$, 3D: $26$), and $\mathbf{H}_{i,i+1}$ is the coupling between the $i$ and $i+1$ shell. We use the superscript $\mathcal{R}$ to indicate the real-space representation of matrices.

We want to calculate the Green function for the infinite matrix comprising the Hamiltonian $\mathbf{H}^{\mathcal{R}}$ in a subspace $\mathbf{H}^{\mathcal{R}}_{0\dots i,0\dots i}\equiv\{\mathbf{H}_{00},\dots,\mathbf{H}_{ii}\}$ up to some shell size $i$. The straight forward Dyson equation is sufficient for systems with short screening lengths such as metals where the convergence requires only a few shellsWu and Jayanthi (1995). For weak screening the increasing matrix sizes with $i$ in the Dyson equation become problematic and one may replace the real-space iterations in shells with an integral over $\mathbf{k}$-points to calculate the real-space Green function (here only shown for the primary unit cell),

[TABLE]

where $\mathbf{G}_{\mathbf{k}}(z)$ is the Green function for a given $\mathbf{k}$-point, $\mathbf{S}^{\mathcal{R}}$ the overlap matrix, and energy plus imaginary part is $z=E+i\eta$. We define the subspace of interest by [math] and the coupling of this to the surrounding bulk system is described by the real-space self-energy, $\boldsymbol{\Sigma}^{\mathcal{R}}_{00}$. We remark that $\mathbf{G}_{\mathbf{k}}^{\mathrm{T}}=\mathbf{G}_{-\mathbf{k}}$ using time-reversal symmetry. This converts the inversion of infinite matrices in real-space to a problem of inverting finite-sized matrices by introducing a $\mathbf{k}$ integral employing Bloch’s theorem. This method was employed in Ref. Li et al., 2018.

Two new problems arise. A sufficient accuracy in the integral is difficult because the elements of the Green function has Lorentzian peaks/step-functions (in $\mathbf{k}$-space) for each eigenvalue (pole). In Fig. 2 we show the Green function matrix elements (left: diagonal, right: off-diagonal) for fixed $z=(0.5+i10^{-4})\,\mathrm{eV}$ using the standard orthogonal tight-binding model for graphene with hopping $t=-2.7\,\mathrm{eV}$. We employ the recursion along one direction in graphene to obtain the Greens function as a function of $\mathbf{k}$ for the direction transverse to this. In effect this means that we sample an extremely dense $\mathbf{k}$-grid along one direction and a sparser $\mathbf{k}$-grid transverse to this. We see that the matrix elements comprise both step-functions and convolutions of Lorentzian and step functions. Such functions require dense integration grids to resolve. Note that using Fourier transforms results in the same deficiencies to resolve the peaks.

A second problem is the matrix dimensions of $\mathbf{G}_{\mathbf{k}}$. Our interest is to calculate the real-space Green function in some multiples of the primary unit cell such that the final matrix has dimensions $n\prod n_{a_{i}}$ with $n$ being the number of orbitals in the primary unit cell, and $n_{a_{i}}$ is the number of repetitions along the $i$’th lattice vector. For large $n_{a_{i}}$ the matrix dimensions rapidly increases making a fine integral in Eq. (3) unfeasibleSettnes et al. (2014).

Our method solves this dimension problem by only performing the costly inversions on matrices of dimension $n$, regardless of all $n_{a_{i}}$. This is achieved using three well established methods and performing them in the following order, 1) a surface self-energy removes the $\mathbf{k}$ integral along a single semi-infinite direction, 2) recursive Green function calculation (block-tri-diagonal inversion, BTD) expands along the semi-infinite direction and 3) Bloch’s theorem efficiently expands the Green function into the remaining one (or two) dimensions. Note that the surface self-energy calculation is a particular efficient solution of the BTD algorithm for a bulk system. Hence the difference between the two methods is subtle, yet significant in their application for our method. The three steps above work for both surfaces and bulk systems, with minor variations. Note that for 3D systems, leaving a $\mathbf{k}$ direction out of the integral, Eq. (3), one finds the real-space Green function for a cylinder with the directions normal to the cylinder surface integrated out, retaining the $\mathbf{k}$ point along the cylinder. Such a particular use-case will not be covered in this study, but we remark that our code allows such calculations which may be useful for e.g. line-defects in solids.

In the following all matrices not denoted by $\mathcal{R}$ are implicitly $\mathbf{k}$ dependent.

II.1 Surface self-energies

The recursive surface self-energy method by Sancho et al. calculates the left/right self-energies for a given transverse $\mathbf{k}$-point. This procedure presents a $2^{i}$ convergence series such that for iteration $i$ one have effectively decimated $2^{i}$ layers. The algorithm is given here for the sake of completeness:

[TABLE]

The surface self-energy removes the $\mathbf{k}$ integral along the semi-infinite direction and immediately reduces the integral in Eq. (3) by one dimension.

II.2 Block-tri-diagonal inversion (BTD)

This method may be generalized to calculate the layer off-diagonals for the inverse of matrices when these can be written in block form Eq. (2) Papior et al. (2017). The pristine bulk system may be written in the following BTD form along the semi-infinite direction:

[TABLE]

Calculating the Green function for an arbitrary number of blocks along the semi-infinite direction follows,

[TABLE]

A key-point is that the real-space Green function for a bulk system is a Toeplitz matrix, e.g. $\mathbf{G}_{mn}=\mathbf{G}_{m^{\prime}n^{\prime}}$ for $m-n=m^{\prime}-n^{\prime}$. Consequently for a bulk system of $M$ blocks one can calculate the full Green function matrix by only calculating $\mathbf{G}_{n0}$ and $\mathbf{G}_{nM}$ for all $n$ (omitting $\mathbf{G}_{MM}$ since it equals $\mathbf{G}_{00}$). Thus only $2M-2$ matrix multiplications are required in order to calculate the full Green function once $\mathbf{G}_{nn}$, $\widetilde{\mathbf{X}}$ and $\widetilde{\mathbf{Y}}$ are obtained. We note that if the system is not bulk (e.g. surfaces) this algorithm need only be replaced by the full BTD algorithmPapior et al. (2017), which is still much faster than full matrix inversion algorithms.

II.3 Bloch’s theorem

We want to obtain the self-energy for the pristine system consisting of a unit cell repeated $N$ times in the transverse direction, large enough to include the defect cf. Fig. 1c. Due to the screening approximation we assume that the potential is unperturbed at the boundary and outside this cell and thus, the self-energy can be calculated from the pristine periodic system. To this end we can apply Bloch’s theorem and express the $N$ times bigger system transverse to the semi-infinite direction via the primary matrix for a given $\mathbf{k}$. In our case we are interested in the Green function for a given $\mathbf{k}$. The equations for expanding the Green function (or any Bloch matrix) along a single direction for a given $K$ (defined in the large $N$ system) is,

[TABLE]

Here $\mathbf{G}_{k_{j}}^{1}$ is the primary cell Green function matrix at the primitive cell $k$-point $k_{j}$ which is to be unfolded into the matrix $\mathbf{G}_{K}^{N}$ and $\otimes$ is the tensor product. The above equation is only expressed in terms of expansion along one direction, however it is easily generalized for more than one direction.

The above three steps conclude the calculation of the real-space Green function for arbitrarily large pristine, periodic systems, $N\times M$. The algorithm in short; the self-energies remove the integral along one $\mathbf{k}$ direction, the BTD algorithm expands the Green function to arbitrary length, $M$, along the semi-infinite direction employing just matrix multiplications, and finally Bloch’s theorem expands the Green function to arbitrary width (and also depth for 3D), $N$.

II.4 Self-energy

While the real-space Green function calculates spectral quantities in a pristine system it is rarely competitive with regular diagonalization methods in the $00$ sub-space and using Bloch’s theorem. Our key mission in calculating the real-space Green function is that it holds the real-space self-energy, $\boldsymbol{\Sigma}^{\mathcal{R}}$, which in turn allows truly single defects (bulk) and contacts (transport) using the Green function formalismBrandbyge et al. (2002); Papior et al. (2017).

In Fig. 3 a schematic calculation shows how the real-space Green function may be used to calculate the real-space self-energy in a region predefined by $M$ and $N$ multiples of the $00$ region as specified in Secs. II.2 and II.3 (hatched region). The real-space self-energy may be conveniently written in two ways:

[TABLE]

where $\mathbf{G}^{\mathcal{R}/00}$ is the real-space Green function for the entire bulk system, excluding the inner region $00$. From Eq. (9) it is clear that $\boldsymbol{\Sigma}^{\mathcal{R}}_{00}$ is non-zero only on sites that connects it to the outside through $\mathbf{V}_{\mathcal{R},00}$. Eq. (10) shows how it is obtained using the real-space Green function.

III Results

In the following we show results on spectral and transport properties of truly single defects/junctions using the real-space self-energy. Our self-consistent DFT$+$NEGF is implemented in Siesta, TranSiesta and TBtransSoler et al. (2002); Papior et al. (2017) while the algorithms described in Sec. II are implemented in sislPapior (2018).

Three systems will be shown using graphene as the real-space electrode. The different systems highlight three particular cases where the real-space self-energy is applicable. We omit the use case of cylindrical self-energies since its use is limited to 3D bulk systems with periodicity along one direction (line defect). The atomic structure of the systems is shown as insets with coloured atoms indicating the support of the real-space self-energy/electrode (in blue), and other electrodes (in red). A last set of atoms is high-lighted (in light green) which are used as the projection region for local density of states (LDOS) analysis.

All calculations are performed using a $300\,\mathrm{eV}$ mesh cut-off, single-$\zeta$ polarized basis set, and PBE$+$GGA exchange-correlationPerdew et al. (1996), and otherwise default parameters. Although LCAO calculations for graphene using simple basis sizes (DZP) misses the lowest unoccupied statesPapior et al. (2018) we do not add basis orbitals to describe these. Thus our presented analysis is limited to energies below the missing unoccupied bands ($E-\mu_{\mathrm{graphene}}<3.35\,\mathrm{eV}$).

III.1 Validation — Graphene

To ensure a functioning method we have constructed a pristine graphene calculation (inset Fig. 4) and calculated the projected DOS on a single carbon atom (marked). The Siesta method calculates the DOS on a $31\times 51$ Monkhorst-Pack gridMonkhorst and Pack (1976) (with energy broadening $\sigma\sqrt{2}=0.1\,\mathrm{eV}$), and both the Green function methods are based on $300\,\mathbf{k}$ points and an imaginary part of $\eta=0.1\,\mathrm{meV}$.

The three different methods all yield the same projected DOS and there is a negligible difference between the methods. Any small difference from the diagonalization method vs. the Green function methods lies in the insufficient $\mathbf{k}$-point sampling and the large smearing value. The two Green function methods are almost numerically the same since the system is bulk and no mixing of $\mathbf{k}$-points take effect.

III.2 Nitrogen defect

Single nitrogen defects in graphene intrinsically have a substantial interaction range and thus calculating defect properties at the DFT level proves difficult Lambin et al. (2012). In this example we use the real-space self-energy and compare with a 2D periodic and 1D periodic calculation. In all 3 examples we use the same unit-cell consisting of a square graphene lattice cell replicated $8\times 9$ totalling 288 atoms. In Fig. 5 we show the projected DOS on the nitrogen atom for the three cases: i) A Siesta calculation using a $31\times 51$ Monkhorst-Pack grid in agreement with other workHou et al. (2013), ii) a two terminal TranSiesta with $300$ transverse $\mathbf{k}$ points, and finally iii) using the real-space self-energy calculated from $300$ $\mathbf{k}$ points. We remark that $300$ $\mathbf{k}$ points corresponds to $2400$ $\mathbf{k}$ in the minimal square graphene unit cell (see e.g. Fig. 2).

The DOS shows distinguished differences and particularly so for energies above the graphene Fermi level. The Siesta and two probe calculations reveals a fine structure with multiple peaks dispersed over $\sim 1\,\mathrm{eV}$. A large smearing parameter ($\sigma\sqrt{2}=0.1\,\mathrm{eV}$) for the Siesta calculation was required due to the relatively crude Monkhorst-Pack grid, which still took more than 5 hours on 20 cores. The two probe calculation shows some even more localized features which could be the same as those in the Siesta calculation. Both look similar to prior calculationsHou et al. (2013) where the projected DOS on the nitrogen defect ($p_{z}$) in a similar periodic simulation was dispersed across two bands with a dispersion $\sim 0.5\,\mathrm{eV}$. We find the real-space method broadens the peaks to a single peak, just above the chemical potential. This result is in perfect agreement with results from a tight-binding description of the isolated N fitted to DFTLambin et al. (2012). Although not shown, the same localized features found for the nitrogen atom are seen for the three neighbouring carbon atoms. These carbon atoms are particularly important for STM imagesLv et al. (2012).

III.3 STM tip on graphene

Scanning tunnelling microscopeBinnig and Rohrer (1983) (STM) is a key experimental technique for analyzing the local electronic structure of surfaces and defects or adsorbates on surfaces. The STM technique is a single tip junction probing the spatial local DOS and yields considerable insight of surface electronic topographies. However, DFT-NEGF calculations of STM on almost isolated defects are problematic both due to periodic repetition of the surface unit-cell, including the repetition of the STM probe tips. Here a calculation of the transmission from an “STM”-like tip to grapheneAndrei et al. (2012); Bhandari et al. (2016); Tetienne et al. (2018); van der Heijden et al. (2016); Márk et al. (2012) is calculated via two methods. Namely, a three terminal (left/right graphene/tip) invoking transverse periodicity, and a two terminal (graphene/tip) calculation, both at an applied bias of $\mu_{\mathrm{graphene}}-\mu_{\mathrm{tip}}=-0.5\,\mathrm{eV}$.

Figure 6 top panel shows spectral DOS of scattering states originating from the tip electrode projected on the carbon atom in contact with tip. In the bottom panel we show the transmission from the tip into graphene.

For the occupied states there is little to no difference while we find a large difference for the unoccupied states. The spectral DOS decreases on the contact atom while the transmission increases. In both graphs we find a discontinuity at $0.45\,\mathrm{eV}$ for the 3-electrode simulation (non existing in the real-space method) which we attribute to periodic image interaction. This fact is supported by other workCalogero et al. (2019) as well as it matches the bias on the tip.

Such STM calculations which may be very influential when calculating inelastic contributionsHalle et al. (2018); Garcia-Lekue et al. (2011) since they use the energies in the unoccupied rangeFrederiksen et al. (2007).

III.4 Graphene contacted to a zGNR

A typical experiment comprise large electrodes contacted through a single junction and rarely arrays of contacts present with few exceptions such as e.g. self-assembled monolayersDubi (2014); Schmaltz et al. (2017). A key issue in DFT$+$NEGF simulations of such devices is that, until now, the simulation had a periodic array of junctions. Such an array of junctions will have interference effects and requires extra care in convergence of the widthThygesen and Jacobsen (2005) and $\mathbf{k}$ points. Using the real-space self-energy we eliminate the periodic junctions and effectively retain a single junction where interference is removed.

The example shown here is a graphene flake contacted to a zig-zag graphene nano-ribbon (zGNR) Leitherer et al. (2019). Our calculations are performed using $\mu_{\mathrm{graphene}}-\mu_{\mathrm{zGNR}}=-0.5\,\mathrm{eV}$. We remark that any molecular junction (for instance Au-benzene-di-thiol-AuMartin et al. (2008); Kergueris et al. (1999)) could be replaced in this example since the electrodes are handled as “surfaces”. In Fig. 7 we plot the projected DOS on the first $4$ atoms in the zGNR (top) and the transmission (bottom). In this example there are relatively few differences since the unit cell is already relatively wide and thus the interference is limited. There are however differences such as a larger spread on the localized states just above the graphene chemical potential. These correspond to states in the zGNR which depends on the electrode coupling and thus is sensitive to periodicitiesPapior et al. (2016); Gunst et al. (2013).

III.5 Electrostatics

For all results shown, the electrostatics has been solved using the Fourier transform. Such a method forces the Poisson solution to be periodic. This is in contrast to our real-space Green function method which inherently has open boundaries. TranSiesta allows external input to ensure the correct boundary conditions for the electrodes. We have ensured that adding such boundary conditions does not change the results noticeably, for further information see (Papior et al., 2017, Sec. 3.5).

III.6 Disk space

Currently, our method relies on storing the self-energies on-disk since the algorithm still needs to be implemented in TranSiesta. Each surface self-energy file has a memory footprint which can be summarized by 3 integer, $n$ (matrix dimension), $n_{\mathbf{k}}$ (number of $\mathbf{k}$-points) and $n_{E}$ (number of energy points)

[TABLE]

where $2+n_{E}$ stems from the Hamiltonian and overlap per $\mathbf{k}$, and the self-energies per $E$. In Tbl. 1 the dimensions and maximum disk usage is shown for the largest bias used.

We find a required disk space requirement of $\sim 50\,\mathrm{GB}$ which is large, but in no way limiting its application on common HPC systems. One generally requires many more energy points in the TBtrans calculation, however since one can define the chemical potential for the real-space electrode to be constant for all applied bias $\mu_{\mathfrak{e}}^{\mathcal{R}}=0$ and the other electrode(s) to be at $\mu_{\mathfrak{e}^{\prime}}=V$ one can reuse the file for all applied bias’ at a much reduced computational cost and with a single file.

III.7 Performance

We have now shown that using the real-space self-energy one can avoid using the super-cell approximation for non-periodic structures such as isolated defects or single junction conductors. In order for it to be competitive with standard methods it also needs to be competitive in terms of performance/through put. We will here show that it is in fact less demanding to do a real-space self-energy calculation when taking into account the full sequence of calculations.

An important factor in using our method is the real-space self-energy calculation. The $\boldsymbol{\Sigma}^{\mathcal{R}}$ method is slower compared to $\boldsymbol{\Sigma}_{L}$/$\boldsymbol{\Sigma}_{R}$ given that the self-energy is more costly to calculate because of larger $\mathbf{k}$-point sampling and a more complex algorithm. On the other hand the SCF cycles and transport/DOS calculations are much faster since no $\mathbf{k}$-point sampling is required. In Tbl. 2 we show the timings of the presented calculations divided into three segments; i) TranSiesta, ii) TBtrans, and iii) $\boldsymbol{\Sigma}^{\mathcal{R}}$. All timings are based on the same 20-core machine.

As can be seen, the timings for TranSiesta is more or less constant while the TBtrans calculations are much faster. Note however, that for the graphene-zGNR system the convergence for the real-space method is faster leading to decreased timings in TranSiesta. Otherwise, the clear bottleneck is the $\boldsymbol{\Sigma}^{\mathcal{R}}$ calculation which can easily be embarrassingly parallelized. We remark that, as noted in Fig. 3, the self-energy is generic for any defect that does not alter the coupling out to the infinite exterior. This means that a single calculation of the self-energy allows using it for more than one system. Since these are one-shot calculations there is no reason not to do an extremely fine $\mathbf{k}$ integration when sampling the real-space self-energy. For the systems shown here it takes less than $100\,\mathrm{s}$ per energy point for $300$ $\mathbf{k}$ points. It should be stressed that the current implementation is done in Python/Cython and thus additional performance gains would be to port it to fortran/C code.

All-in-all we find that the proposed method is comparable to, or faster, than the existing method for equivalent $\mathbf{k}$-point sampling.

IV Conclusion

We have presented a simple, effective, and fast algorithm for constructing real-space self-energies generalized for surfaces and full 2D/3D bulk systems. The algorithm relies on already well established methods used in the community and can thus be directly integrated into existing codes without problems. The current algorithms are implemented in the TranSiesta, TBtrans and sisl toolboxes which are all open-source under GPL variant licenses.

We have applied the method in three graphene cases which are readily found in current experimental literatureLv et al. (2012); Andrei et al. (2012); Halle et al. (2018); Martin et al. (2008); Kergueris et al. (1999). A recurring difference between the analyzed DOS and transmission profiles is that the occupied energy range is largely comparable to standard DFT$+$NEGF methods, while the unoccupied energy range shows substantial deviations. Such differences are attributed to removed interference effects.

We have shown how the use of real-space self-energies will remove the periodic images of defects in DFT calculations. The results shown provide insights into the far-field accuracy of DFT$+$NEGF calculations for single defects which has been missing in the electronic structure community.

V Acknowledgments

We thank Dr. K. Kaasbjerg for useful discussions. Financial support by Villum Fonden (00013340), Danish research council (4184-0003). The Center for Nanostructured Graphene (CNG) is sponsored by the Danish Research Foundation (DNRF103).