Electronic Transport (NEGF)¶

Some examples are available in the ADFHOME/examples/band directory and are discussed in the Examples section. Note In the BAND-GUI it is possible to choose between three NEGF methods (flavors): Self consistent This is the internal BAND-NEGF implementation, which is described in this page. Self consistent + align This is the internal BAND-NEGF implementation with an extra alignment-run (workflow step 3a) Non self consistent Computationally cheap method, equivalent to the DFTB-NEGF approach with $$H$$ and $$S$$ matrix elements computed by BAND (instead of DFTB). Transport with NEGF in a nutshell¶ The Non-Equilibrium Green’s Functions formalism (NEGF) is a theoretical framework for modeling electron transport through nano-scale devices. Electron transport is treated as a one-dimensional coherent scattering process in the “scattering region” for electrons coming in from the electrodes: Our goal is to compute the transmission function $$T(E)$$, which describes the rate at which electrons of energy $$E$$ are transferred from the left electrode to the right electrode by propagating through the scattering region. From the transmission function we can calculate the electric current for given Bias Voltage $$V$$ applied between the electrodes: $I(V) = \frac{2e}{h} \int_{-\infty}^\infty T(E,V) \left( f(E - \mu_L) - f(E - \mu_R) \right) dE$ where $$f(E)$$ is the Fermi-Dirac distribution function for a given temperature, and $$\mu_L$$ ($$\mu_R$$) is $$\epsilon_F + eV/2$$ ($$\epsilon_F - eV/2$$), $$\epsilon_F$$ being the Fermi energy of the electrodes. The transmission function $$T(E)$$ can be computed from the Green’s function of our system. The Green’s function $$G(E)$$ of the scattering region is obtained solving the following equation: $(E S - H) G(E) = I$ where $$S$$ is the overlap matrix, $$H$$ is the Hamiltonian and $$I$$ is the identity matrix. The Hamiltonian is composed as follows (L, C and R denote the left lead, the central region and the right lead respectively): $\begin{split}H = \left( \begin{array}{ccc} H_L + \Sigma_L & H_{LC} & 0 \\ H_{LC} & H_C & H_{RC} \\ 0 & H_{RC} & H_R + \Sigma_R \end{array} \right)\end{split}$ The two self-energies $$\Sigma_L$$ and $$\Sigma_R$$ model the two semi-infinite electrodes. The transmission function $$T(E)$$ can be calculated from the Green’s function $$G(E)$$ and the so-called broadening matrices $$\Gamma_L(E)$$ and $$\Gamma_R(E)$$: $T(E) = Tr[G(E) \Gamma_R(E) G(E) \Gamma_L(E)]$ The broadening matrix being $\Gamma_L(E) = -2 \Im \Sigma_L(E)$ Self consistency¶ The density matrix is determined self consistently [68]: P_\text{in} \rightarrow H_{KS} \stackrel{\mbox{shifts}}{\longrightarrow} H_\text{aligned} + \Sigma_L(E) + \Sigma_R(E) \rightarrow G(E) \stackrel{\int de}{\longrightarrow} P_\text{out} From a guess of the density matrix the corresponding KS Hamiltonian is calculated. This Hamiltonian is aligned, and then the NEGF Hamiltonian in the complex plane is constructed by adding the self energies, representing the influence of the electrodes. From the resulting Green’s function a new density matrix follows. From the difference between input and output density a next input is guessed. This is repeated until the input and output densities converge. For the alignment of the Hamiltonian there are two shifts. The first shift aligns the potential in the leads to the electrodes. $\mbox{shift 1} = \frac{1}{n} \sum_\text{i in lead}^n \frac{H^{TB}_{ii} - H^{KS}_{ii}}{S_{ii}}$ The second and usually smaller shift results from the alignment run. A shift $$\Delta$$ is applied globally H_{ij}^\text{aligned} = H_{ij} + \Delta S_{ij} Contour integral¶ Without bias the density matrix follows from $P(\mu) = -\frac{1}{\pi} \int_{-\infty}^{\infty} de \; f(e,\mu) \; \Im G(e)$ As the Green’s function is singular on the real axis we add a small imaginary value (eta) to the energy. Still, the integrand will be very wild function, and it is numerically better to do a contour integral instead. Gate potential¶ There is no direct key for the gate potential. You can model this with the FuzzyPotential key. Setting up the gate potential for NEGF is most conveniently done with the GUI. Bias potential¶ When there is a bias specified there are two important things to keep in mind. First of all you need to define a ramp potential. In the negative lead this should have the value +V/2 and in the negative lead -V/2. The ramp should smoothly go from one to the other value. For metals one could start the ramp at the surface atoms of the lead material. For semi-conductors it is less clear. The ramp potential can be specified with the FuzzyPotential key. The GUI can be helpful here. Secondly, the expression for the density is different from the zero-bias case: $\rho = \rho_\text{eq}(\mu_+) + \rho_\text{neq}$ The first (equilibrium) term is calculated with a contour integral as before, the second (non-equilibrium) part cannot be calculated with a contour integral. Instead, an integral in the complex plane (close to the real axis) is performed, the range covering the bias energy window. See also PhD Thesis of C. Verzijl (BAND-NEGF developer) Workflow¶ The computation of the transmission function $$T(E)$$ within the BAND-NEGF [68] formalisms requires three or four individual simulations. Tip Use ADFInput (GUI) to set up your BAND-NEGF calculation (see the BAND-NEGF GUI tutorial) 1): Lead calculation A 1D-periodic BAND calculation of the lead (including StoreHamiltonian2): A tight binding (TB) representation is calculated for the overlap ($$S(R=0)$$ and $$S(R=a)$$) and Fock matrix ($$H(R=0)$$ and $$H(R=a)$$). This is not an approximation provided that the functions do not extend beyond the neighboring cells. You should choose a sufficiently large super cell for this to be true. For this reason we recommend setting the SoftConfinement Quality to Basic, thus reducing the range of the functions. 2): SGF calculation A small program that determines the fermi energy $$\epsilon_F$$ corresponding to the TB representation, and the specified temperature. This fermi energy is typically a bit higher than the one from the lead calculation. This also tests the contour integration. 3a): Alignment run (optional) The idea is to fill the central region with bulk material. Then one expects to have zero charge in the central region. In practice this is not exactly true. In the alignment run the shift is determined that makes the central region neutral. This global shift is to be used in the next run. 3b): Transport calculation Computes the NEGF transmission function $$T(E)$$. The density matrix is determined fully self-consistently. Without alignment (3a) one should set NEGF%ApplyShift2 to False. To get the current as a function of bias potential you need to repeat calculation 3b for a various bias potentials. SGF Input options¶ SGF is a small separate program. An input looks like: ADFBIN/sgf   << eor
TITLE Test for NEGF inputs
SAVE SIGMA
SURFACEGF
SCMCode
KT 0.001
ContourQuality normal
END
eor


It looks for a file RUNKF and the output is a file named SigmaSCM. The only important parameter is KT which is the Boltzmann constant times the temperature in Hartree. The other parameter of interest is the ContourQuality, which can be set to Basic,Normal,Good,VeryGood, or Excellent.

NEGF Input options (no bias)¶

The NEGF functionality is controlled by the NEGF block key:

NEGF
SGFFile Sigma.runkf
ContourQuality [Basic|Normal|Good|VeryGood]
EMin -0.1837465475
EMax 0.1837465475
NE 200
! From here one the options are very technical
CheckOverlapTol rval
Eta rval
ApplyShift1 True
ApplyShift2 False
yContourInt rval
deContourInt rval
End

LeadFile
Contains the tight binding representation of the lead.
SGFFile
The result from the SGF program. Contains the fermi energy of the lead.
ContourQuality
Default: normal. The density matrix is calculated numerically via a contour integral. Setting this to Basic,Normal,Good,VeryGood influences the number of points. This influences a lot the performance.
NE
The nr. of energies for the transmission energy grid.
EMin
The minimum energy for the transmission grid (with respect to the fermi level of the lead).
EMax
The maximum energy for the transmission grid.

The remaining options are very technical in nature.

CheckOverlapTol
Default: 0.01 BAND checks how well the TB overlap matrix :math:S(R=0) represents the overlap matrix in the lead region. Elements corresponding to the outer layer are neglected, because when using a frozen core they have bigger errors.
ApplyShift1
Default: True Apply the main shift, obtained from comparing matrix elements in the leads with those from the tight-binding run. Strongly recomended.
ApplyShift2
Default: True Apply the smaller alignment shift. This requires an extra alignment run. Usually this shift is smaller.
Eta
Default: 1e-5 Small value used for the contour integral: stay at least this much above the real axis. This value is also used for the evaluation of the Transmission and dos.
yContourInt
Defaults: 0.3 The density is calculated via a contour integral. This value specifies how far above the real axis the (horizontal part of the) contour runs. The value is rounded in such a way that it goes exactly halfway between two fermi poles. There is a trade off: making it bigger makes the integrand more smooth, but the number of enclosed poles increases. For low temperatures it makes sense to lower this value, and use a smaller deContourInt.
deContourInt
The energy interval for the contour grid. Defaults depends on the contour quality.

NEGF Input options (with bias)¶

With a bias potential there are some extra keys.

NEGF
BiasPotential rval
NonEqDensityMethod [1|2|3]
BoundOccupationMethod [1|2|3]
yRealAxisInt rval
deRealAxisInt rval
End

BiasPotential
Apply a bias potential (atomic units). Can be negative. One has to specify the ramp potential with the FuzzyPotential key. This is mostly conveniently done with the GUI.
NonEqDensityMethod

Default: 1 Let us introduce some terms [66]. First of all the total density in the bias window (ignoring occupation)

$D = \frac{1}{2\pi} \int A \; (f_--f_+)$

And then there are the side resolved densities

$D_{+/-} = \frac{1}{2\pi} \int A_{+/-} \; (f_--f_+)$

The issue here is that the side resolved densities do not sum to the total one

$D = D_{+} + D_{-} + D_\text{bound states}$

The NonEqDensityMethod is about how these integrals are calculated . With option 1, or 2 a contour integral is used for D: they are essentially the same. However, when choosing option 2, you can choose a BoundOccupationMethod, leading to other physics. If set to 3, the total density in the bias window (D) will be calculated near the real axis: this way one avoids the possibility of a negative nr. of bound states (deviating from [66]).

BoundOccupationMethod

Default: 1 Only relevant with NonEqDensityMethod equal 2 or 3. If set to one, the density of bound states (ignoring occupation) is simply multiplied by a half. If set to two, atoms closer to the negative lead will get a higher occupation [66]. Atoms coupled to the right lead will have a low occupation. For this we recommend setting NonEqDensityMethod to 3, to avoid a possible negative number of bound states.

Setting the method BoundOccupationMethod to 1, leads to

$\rho = \rho(\mu_+) + \rho_\text{-/neq} + \frac{1}{2}\rho_\text{bound}$

By setting the method to 2, each atom gets its own weight in the density matrix

$\rho_{ij} = \rho_{ij}(\mu_+) + \rho_\text{-/neq} + \sqrt{w_i w_j} \rho^\text{bound}_{ij}$

with [66]

$w_i = \frac{\text{Tr} [D_-]_i}{\text{Tr} [D_-]_i + \text{Tr} [D_+]_i}$

These weights are the same for all functions on an atom. The intended effect is: bound states that are coupled more strongly to the negative electrode get a higher occupation than the ones that are coupled more strongly to the positive electrode.

To summarize here are three reasonable settings

 NonEqDensityMethod BoundOccupationMethod intention 1 1 Multiply the bound states with a half 2 2 Occupy bound states with atom-resolved $$w_i$$ 3 2 ... and prevent a negative nr. of bound states
yRealAxisInt
Defaults: 1e-5 The non-Equilibrium density is calculated near the real axis.
deRealAxisInt
The energy interval for the real axis grid. Defaults depends on the contour quality.

To get the current from a calculation you can use adfprep:

\$ADFBIN/adfreport RUNKF 'NEGF%current'


NEGF Input options (alignment)¶

For the (optional) alignment run there are some extra keys.

NEGF
DoAlignment True
Alpha rval
AlignChargeTol rval
CDIIS [True|False]
End

DoAlignment
Default: False. Set this to True if you want to do an align run. Between the leads there should be lead material. The GUI can be of help here.
Alpha
Default: 1e-5. A charge error needs to be translated in a potential shift.
$\Delta V = \alpha \Delta Q$
AlignChargeTol
Default: 0.1 In an alignment run you want to get the nr. of electrons in the center right. This number specifies the criterion for that.
CDIIS
Default: False. Make the normal DIIS procedure aware of the align charge error.

Troubleshooting¶

The self consistent approach, unique to BAND, may be difficult to converge. If this is true for the alignment run it can be decided to skip this run. For the final transport run, here are some tips / considerations.

• Use a SZ basis for the metal atoms
• Restart (the density matrix) from the result of a smaller (such as the SZ) basis. (See “Save DensityMatrix” and the Restart key)
• Restart (the density matrix) from the result obtained with a smaller bias (only relevant for calculations with bias potential).
• Setting NEGF%BoundOccupationMethod to 2 (and NEGF%NonEqDensityMethod to 3) might help. Note that this affects the physics: you are differently occupying the bound states.
• Use a better NEGF%ContourQuality (there comes a computational price tag with this).

If everything fails it is possible to use BAND in a non-self consistent way, similar to the way DFTB-NEGF works. This option is available via the GUI.

Miscellaneous remarks on BAND-NEGF¶

• You should make sure that your results are converged with respect to the number of lead repetitions; the results should not change significantly if you increase the number of lead repetitions.
• It’s good practice to include at least one lead repetition in the central region.

Store tight-binding Hamiltonian¶

Let us consider a Fourier transformation of a 1D Bloch matrix

$S(R=na) = \int_k e^{-ikR} S(k)$

In (the tight-binding) case that the functions do not extend beyond the neighboring cells only S(R=0) and S(R=a) are nonzero. (And S(R=-a) is equivalent to S(R=a))

StoreHamiltonian2


Adding StoreHamiltonian2 to the input cause band to determine the tight-binding representation of the overlap an fock matrix. Currently this only works for 1D periodic systems. For the overlap matrix you will get two parts. The first $$S(R=0)$$ is the (symmetric) overlap matrix of atoms in the unit cell. The second $$S(R=a)$$ is a non symmetric matrix describing the coupling of functions in the central cell with functions in its right neighboring cell. On the RUNKF file you will find the TB representations of the overlap and Hamiltonian stored in the ‘Matrices’ section as “S(R)” and “H(R)”, being dimensioned (nBas,nBas,2).