Non-equilibrium Green's functions in density functional tight binding: method and applications

We present a detailed description of the implementation of the non-equilibrium Green's function (NEGF) technique on the density-functional-based tight-binding (gDFTB) simulation tool. This approach can be used to compute electronic transport in organic and inorganic molecular-scale devices. The DFTB tight-binding formulation gives an efficient computational tool that is able to handle a large number of atoms. NEGFs are used to compute the electronic density self-consistently with the open-boundary conditions naturally encountered in quantum transport problems and the boundary conditions imposed by the potentials at the contacts. The efficient block-iterative algorithm used to compute the Green's functions is illustrated. The Hartree potential of the density-functional Hamiltonian is obtained by solving the three-dimensional Poisson equation. A scheme to treat geometrically complex boundary conditions is discussed, including the possibility of including multiterminal calculations.

3 function computationally efficient. In section 5, we describe the Poisson solver used to calculate the Hartree potential in the molecular region, needed to compute the non-equilibrium charge density. Applications to molecular devices in three-terminal FET configurations are shown in sections 6 and 7.

The DFTB method
The DFTB approach consists of constructing a non-orthogonal tight-binding (TB)-like representation of the single-particle wavefunctions, starting from a DFT formulation of the many-electron problem. The atomic wavefunctions are calculated from one-center all-electron DFT calculations with appropriate compression potentials [13].
Like in empirical TB, the system Hamiltonian is not calculated on the fly, but is obtained from numerical look-up tables. The Hamiltonian and overlap matrix elements between pairs of atomic wavefunctions are numerically pre-calculated using ab initio DFT over a grid of inter-atomic distances. Three-center integrals are neglected [11]. The method is also extended to self-consistent K-S calculations using an approximate local density approximation (LDA) exchange-correlation functional [16]. In this scheme, the charge density is expanded on a reference density and a correction, n = n 0 + δn, and the K-S equations result in where H 0 is the TB Hamiltonian obtained for the reference density n 0 , H 1 is the self-consistent correction of order O(δn), S is the overlap matrix between orbitals, and C i and E i are, respectively, the eigenvectors and eigenvalues of the secular equation.
The self-consistent part can be expressed as [16] H 1 µν = 1 2 S µν k (γ ik + γ jk ) q k , ∀µ ∈ i, ν ∈ j, (2) where q k are the Mulliken charges and γ ik are appropriate functionals taking into account Hartree and exchange-correlation interactions, where the U H i are the Hubbard on-site repulsion energies and the F i 00 (r ) are appropriate spherical densities centered around the atoms [16]. Since the atomic charges depend on the oneparticle wavefunctions, a self-consistent procedure is required. The improvement of the selfconsistent over the non-self-consistent procedure is considerable in determining the structural and energetic properties of molecular systems [16].
The TB formulation of the electronic energy is complemented by the semi-empirical construction of inter-atomic repulsive potentials, which allows the computation of atomic forces and enables molecular dynamics (MD) calculations of large super-cells. This is particularly suitable for studing the electronic properties and dynamics of large systems such as proteins, carbon nanotubes, DNA strands or adsorbates on surfaces, semiconducting heterostructure, etc (see [17] for applications).

The NEGF extension
Despite its mathematical complexity, the NEGF method for calculation of quantum transport has gained great popularity in recent years, mostly because of the versatility and numerical stability of the method, in contrast to wavefunction or transfer matrix approaches. The open boundary conditions can be elegantly included by exactly mapping the contacting leads into a finite and small part of the system [18]; furthermore, the Green's function approach can be generalized to many-body quantum theory, allowing the inclusion of electron-phonon [19,20] as well as electron-electron interactions [21]- [23] within a unified and systematic formalism. Good references on many-body quantum theory can be found in [24,25] and an exhaustive review on NEGF can be found in [26].
The types of systems under study can be represented as in the graph of figure 1. The system can be arbitrarily partitioned into the contacts C α and a device region, D. An arbitrary number of contacts can be included. The contacts are semi-infinite leads and it is assumed that their properties coincide with those of bulk systems [27]. The device is a collection of atoms linking the contacts, comprising a molecular bridge, M, and a portion of the surfaces of the conducting leads, S α . The inclusion of the surfaces within the molecular region is necessary in order to ensure that the regions C α can really be considered bulk-like. This assumption can be directly verified by checking that the charge density smoothly joins at the boundaries C α /S α . The molecule plus the surfaces can be called extended molecule. The contacts are kept at different electrochemical potentials, driving current fluxes across the device.
In solving for the transport problem, it is not possible to make the assumption of local thermodynamical equilibrium, i.e. a global Fermi energy is not defined. The only assumption that can be made is that the connecting leads are kept at different electrochemical potentials, and are considered large reservoirs where the electrons are effectively in equilibrium. In order to be consistent with this scenario, the conducting bridge must offer the largest source of resistance to the flowing current. Only under this condition is it consistent to assume that the potential drops essentially across the device region, while the contacting leads are in equilibrium at two different constant potentials. The density matrix is computed, as usual, as where G < is the electron-electron correlation matrix, represented on the local basis, and is proportional to the spectral density of occupied states. Another important quantity is G > (E), representing the density of empty states. The spectral density of states can be expressed in terms of the retarded Green's function as It is important to remark that all results derived here are valid in steady-state conditions, such that the two-time Green's functions, G(t, t ), only depend on the time difference, t − t , and can be Fourier-transformed to energy. In equation (6), we have explicitly included the retarded self-energy, r (E), which accounts for the presence of the contacts and can account for other perturbing potentials due to phonons, impurities or the other electrons.
One of the advantages of this approach is that the contacting leads can be exactly mapped into the extended molecule using appropriate self-energies (at least within the single-particle approximation). In local orbital representations, where the interaction between atoms has a finite range, the contact self-energy can be easily calculated by exploiting the fact that the Hamiltonian describing the interaction between the device region and the contacts involves a finite number of atoms close to the junctions. Therefore, the required contact Green's function can be solved just for the matrix block corresponding to the atoms close to the extended molecule region. This so-called surface Green's function can be calculated using powerful recursive algorithms [28]. The density matrix can be used to compute the Mulliken charges, according to where q 0 i are the reference atomic charges usually chosen as neutral atom charges. The dynamics of the occupation due to inelastic scattering processes from one energy channel to the other is provided by the Keldysh-Kadanoff-Baym (KKB) equation [29,30], also known as kinetic equation. Within the TB matrix representation, the steady-state solution of the KKB equation can be expressed in terms of matrices, as [31] where G a = G r † is the advanced Green's function, and < and > are the non-equilibrium self-energies.
In practical computations, it is much more convenient to divide G < into three contributions in order to obtain the density matrix, where f c is related to the collector contact, i.e. the contact at the lowest electrochemical potential. The first term can be efficiently computed on a complex contour [32]. The contour is deformed in the upper complex plane since the poles of G r are slightly displaced below the real axis. The second term is obtained from the result of the first by taking the adjoint. Finally, the integration involved in the third term must be performed along the real axis, since the integrand is analytic on the real axis only. Here we have defined, as usual, . Although this integration is restricted by the Fermi functions to a segment of the order of the applied bias voltage, it may be delicate and requires a fine mesh.

The block-iterative algorithm
In local basis representations, the Hamiltonian and the overlap matrices have a sparse structure because of the short inter-atomic interactions. In quasi-one-dimensional systems, the atoms can be arranged in layers, such that H and S take the block tridiagonal form represented in figure 2. The diagonal blocks can be of different sizes, but all of them must have a square shape. The corresponding layers are called principal layers (PLs) and these have the property that the interaction of each PL is restricted to two nearest-neighbor PLs. This structure can be exploited by devising an efficient computation of the Green's function. Indeed, the direct inversion of the whole Hamiltonian can be extremely time-consuming in large structures. The equation defining the Green's function, can be written in terms of a tridiagonal Hamiltonian. We conveniently define the off-diagonal terms as 7 In general, the Green's function does not preserve the block tridiagonal structure, but only the sub-blocks where the matrix S is non-vanishing are actually needed to compute the Mulliken charges, thanks to equation (7). The iterative implementation described here has been optimized for one-dimensional topologies. Imagine we deal with a device comprising N layers, interacting with two or more contacts. We assume that two such contacts interact with the device at layer 1 and at layer N . The remaining contacts may interact with any of the intermediate layers. As a preliminary step, the contact self-energies are computed, according to the equation where L(α) is the layer of the device interacting with contact α, and g α,α is the surface Green's function of contact α. The latter is computed assuming that the contact Hamiltonian can also be represented in a block tridiagonal form, with all diagonal blocks and interactions being the same. are obtained by removing the interactions T L ,L+1 and T L+1,L . In practice, g L L ,L is the exact surface Green's functions of a system cut at layer L and interacting only with the left part of the system. On the other hand, g R L ,L is the surface Green's functions of a system cut at layer L and interacting only with the right part of the system. These Green's functions obey the following recursive relations: and The Green's functions g R L ,L are computed by iterating upward, starting from layer N up to layer 1, whereas the Green's functions g L L ,L are computed by iterating downward, starting from layer 1.
In our implementation, we compute first the partial Green's functions, g R L ,L , then we can start computing the complete Green's functions iterating downward, starting from layer 1, and then all subsequent layers down to N , using The off-diagonal blocks can be obtained from the expressions and All the expressions derived so far are valid for G r or G a . In equilibrium, these are sufficient to compute the density matrix. In order to derive the correct expressions for the non-equilibrium case, we need to carry the matrix multiplication (8) explicitly. Again, we only need the tridiagonal blocks of the matrices G <,α = G r α G a , which can be written as This set of equations shows that the computation of G <,α L ,L requires the computation of the column-blocks G r L ,L(α) , corresponding to the layer L(α). Referring back to equation (9), we observe that for a system of N contacts, the computation of N − 1 columns is needed. These column-blocks can be obtained using the recursive formulae valid for j > L(α), and valid for j < L(α).
A further complication in the computation of the Mulliken charges is the non-zero overlap between device and contact atoms. As a consequence of equation (7), the density matrix needs to be computed on extra blocks outside the device, corresponding to the overlap between device and contacts, as shown in figure 2. Concerning G r , the calculation of these blocks is simple and can be derived, in analogy to equation (17), as Things are a little more complicated for the corresponding blocks of G <,α . With lengthy algebra, it is possible to obtain the general expressions and  Figure 3 shows the comparisons between the computational time and memory load of a direct inversion of the whole matrix and the block iterative algorithm applied to a Si nanowire comprising identical layers of 60 Si atoms each, corresponding to a block matrix 240×240 using an sp 3 basis set. It can be seen that the iterative algorithm scales linearly with the wire length both in computational time and peak memory load.

The Poisson equation
As anticipated in section 3, the Hartree potential needed for the self-consistent iteration of the K-S equations is computed by solving the Poisson equation with appropriate boundary conditions imposed by the contact potentials. The total electronic density, n(r ), can be split into a reference density of neutral atoms and a density fluctuation as where n 0 i (r ) are the atomic reference densities, F i 00 (r ) are spherical s-like radial functions and q i are the Mulliken charges. The Poisson equation for the mean-field electrostatic potential should be written as In DFTB, the self-consistent potential is related to the electronic density fluctuations, whereas the effective potential, V 0 , corresponding to the reference charge density, n 0 (r ), is included in H 0 . Hence, by linearity, we split the Poisson equation (28) into two equations, one for the reference density and the other for the self-consistent correction, When solved with the usual boundary conditions (potential vanishing at infinity), equation (29) gives the usual Hartree potential, This potential can be projected on the atom centers by using which can be used to recover the gamma functional of equation (4).
In gDFTB, equation (29) is solved with the boundary conditions imposed by the contact potentials. These conditions arise from the natural requirement that deep inside the contacts the effective potential for the K-S equations must correspond to the bulk electrochemical potentials. Therefore, at the boundaries between the device region and the contacts, the potential must match the intrinsic effective bulk potential (which originates from any equilibrium charge density) shifted by the applied bias. At the device-contact interfaces, C α /S α , the potential must satisfy where V α is the applied external potential to the α-contact. The decoupling of equation (28), which at first may seem arbitrary, is actually a good approximation since the reference density is taken as that of the neutral atoms and therefore cancels with the ionic charges. In contrast, the excess density produces a long-range Coulomb field that should respect the boundary conditions imposed by the device. For example, the charge which accumulates on the contact surfaces must be consistent with the applied bias.
Within the gDFTB approach, the Poisson equation is solved in real space using a threedimensional multigrid algorithm applied to a general linear, non-separable, elliptical partial differential equation. The Poisson equation is just a particular case of this kind of equation, which has the form i=x,y,z The The mesh is usually chosen as a trade-off between accuracy and computational speed. The usual spherically symmetric form used in DFTB for the atomic charge density [16] gives quite smooth functions, and usually convergent results are obtained with a mesh spacing of 0.5 atomic units.

A donor-acceptor rectifying molecule
In this section, we consider a donor-acceptor-type molecule, represented in figure 4, comprising a tetrathiafulvalene (TTF) unit, a benzo-quinone-diimine (BQD) unit, separated by a phenyl ring. An ethyl-dioxy-thiophene (EDT) group has been attached to the central benzene ring for reasons that will be clear later. The molecule is first relaxed in a free-standing configuration, with S atoms saturated with H. Subsequently, the terminating hydrogens are removed and two Au contact leads are placed at a distance of 2.5 Å from the S atoms, which is known to be the equilibrium distance from previous relaxations of thiols on Au surfaces. The Au leads are taken as Au(110) nanowires that mimic two tips.
The TTF acts as an electron donor with the HOMO energy level at −4.2 eV and the LUMO level at −2.1 eV. The BQD side acts as an electron acceptor with the HOMO level at −4.8 eV and the LUMO level at −2.8 eV. The benzene ring in between acts as an insulating bridge that separates the two molecular islands, because of its larger energy gap. In order to increase localization of the states on the two dots, we have introduced a CH 2 group between the phenyl ring and the BQD. We define forward bias as a positive potential applied to the contact at the BQD side, while the TTF side is grounded. In forward bias, the HOMO level at the TTF side is brought into resonance with the LUMO level of the BQD side, leading to a state of high The potential is drawn on a plane cutting the molecule and is referenced with respect to V = 0. Black atoms correspond to C, green atoms to S, blue atoms to N and yellow atoms to Au. conduction of the molecular bridge, as shown in figure 5. In reverse bias, the same levels are moved far apart, giving a low tunneling current. This rectifying behavior is shown in the I -V characteristics plotted in figure 6. Figure 4 shows the potential drop between the two contacts at the bias of 1.5 V. This potential is obtained as a difference with respect to the reference potential at V = 0 in order to eliminate local atomic fluctuations. The plot confirms that the largest potential drop occurs across the benzene ring, as expected.  The molecular structure is shown in the inset.

A three-terminal FET-like configuration
In this section and the next section, we present the first application of a non-equilibrium calculation in a three-terminal configuration in which the potential is solved fully selfconsistently and the tunneling current flowing between arbitrary pairs of contacts can be analyzed in detail. The TTF-EDT-BQD molecule, studied in section 6, has been arranged in a three-terminal configuration, as represented in figure 7. A sulfur atom has been added to the EDT group to bind to a third Au-wire. The molecule is finally relaxed in the three-terminal configuration, keeping the Au atoms fixed. The three-terminal configuration can be thought of as a molecular analog of a FET transistor, where the TTF-BQD acts as a source-drain channel and the contact attached to the EDT group plays the role of a gate that modulates the conductance. The gate field modulates the HOMO-LUMO energy alignment of the central benzene ring with respect to the FFT and BQD levels, increasing or decreasing the tunneling current. By applying a negative potential to the gate, the HOMO level of the benzene ring is moved upwards, getting closer to the HOMO-LUMO levels of the TTF-BQD and favoring the tunneling current. The drain bias is applied to the BQD side of the molecule. Figure 8 shows the I -V characteristics of the three-terminal device, at the gate voltages of 0, −0.2, −0.5 and −0.7 V. We observe a certain degree of modulation induced by the gate field, however, for drain biases larger than 0.2 V, the drain current starts to decrease considerably, and the modulation effect is considerably reduced for drain voltages larger than 1.0 V. This result can be understood in terms of the socalled 'drain-induced barrier lowering' known in semiconductor device physics. In practice, the gate field is not able to pin the potential of the benzene ring, which is strongly affected by the drain voltage. The effect can be clearly seen in figure 9, where the peak corresponding to the benzene ring is strongly shifted downwards by the drain potential rather than being pinned by the gate. Another drawback of this molecular configuration is the rather high source-gate leakage current, which overcomes the source-drain current by up to one order of magnitude. This is due to the good degree of conjugation between TFF and EDT, which is broken in the TTF-BQD by the alkyl group. This problem can be overcome for instance by increasing the degree of conjugation of the 'channel', e.g. by replacing the alkyl separation between the central phenyl ring and the BQD with a sulfur atom, and by placing an alkyl group to improve the isolation of the EDT island from the channel, as shown in figure 10. This molecular structure has been obtained with the same procedure followed for the previous structure. The result of this strategy is plotted in the  I -V curves of figure 11, showing a much larger output current with respect to figure 8 owing to a larger degree of conjugation and a slightly better gate pinning, particularly at the highest gate voltages. Also the gate leakage, as shown in figure 12, is significantly reduced in this molecular configuration thanks to the alkyl group. The reduction of the leakage current is indeed about one order of magnitude.

Conclusions
In this paper, we have reviewed the NEGF theory applied to the DFTB code for the calculation of self-consistent potentials in nano-and molecular junctions. We have discussed the implementation of an iterative Green's function algorithm and a multigrid Poisson solver used to efficiently calculate complex geometries, comprising several hundreds of atoms and many contacts. The resulting code, known as gDFTB, has been applied to several molecular systems in two-and three-terminal configurations. In particular, we have analyzed the I -V characteristics of a tetrathiafulvalene-phenyl-benzoquinone diimine (TTF-EDT-BQD) rectifier and a TTF-EDT-BQD in a three-terminal device used as a field-effect transistor. From our calculations we learn that FET operation depends on several parameters, including the topological configuration of the molecular subunits with respect to the contacts, which affects the gate and drain control. Furthermore, we see that it is not easy to obtain good gate pinning, because the drain field can easily dominate the molecular electrostatics. A good compromise between high molecular conduction and isolation of the gate is necessary for good device operation.