Efficient implementation of the nonequilibrium Green function method for electronic transport calculations

An efficient implementation of the nonequilibrium Green function method combined with the density-functional theory, using localized pseudoatomic orbitals, is presented for electronic transport calculations of a system connected with two leads under a finite bias voltage. In the implementation, accurate and efficient methods are developed especially for the evaluation of the density matrix and treatment of boundaries between the scattering region and the leads. Equilibrium and nonequilibrium contributions in the density matrix are evaluated with very high precision by a contour integration with a continued fraction representation of the Fermi-Dirac function and by a simple quadrature on the real axis with a small imaginary part, respectively. The Hartree potential is computed efficiently by a combination of the two-dimensional fast Fourier transform and a finite difference method, and the charge density near the boundaries is constructed with a careful treatment to avoid the spurious scattering at the boundaries. The efficiency of the implementation is demonstrated by rapid convergence properties of the density matrix. In addition, as an illustration, our method is applied for zigzag graphene nanoribbons, a Fe/MgO/Fe tunneling junction, and a ${\text{LaMnO}}_3/{\text{SrMnO}}_3$ superlattice, demonstrating its applicability to a wide variety of systems.

Figure 1

(a) Configuration of the system, treated by the NEGF method, with infinite left and right leads along the $\mathbf{a}$ axis under a two-dimensional periodic boundary condition on the $\mathbf{b}\mathbf{c}$ plane. (b) One-dimensional system compacted from the configuration of (a) by considering the periodicity on the $\mathbf{b}\mathbf{c}$ plane, where the region $C$ is an extended central region consisting of $C_0$, $L_0$, and $R_0$.

Figure 2

(Color online) (a) Absolute error in $\left(E_{\text{tot}}-E_{\text{self}}\right)$ per atom of a linear carbon chain with a bond length of $1.4\text{Å}$ under the zero-bias voltage at electronic temperature of 300, 600, and 1200 K, where the regions $L_0$, $R_0$, and $C_0$ contain four carbon atoms, respectively. For comparison, the same calculation (closed contour) using a closed contour method (Refs. 16, 67) is also shown for $T=600\text{K}$. The definitions of $E_{\text{tot}}$ and $E_{\text{self}}$ are found in Appendix . The reference values are obtained from calculations with a large number of poles. (b) Total DOS of the carbon linear chain, calculated by the conventional band-structure calculation (solid line) and the EGF method (dotted line), under the zero-bias voltage at 300 K, where 160 poles are used for the integration of the equilibrium density matrix. It is hard to distinguish two lines due to the nearly equivalent results.

Figure 3

(Color online) (a) Absolute error in $\left(E_{\text{tot}}-E_{\text{self}}\right)$ per atom of the same linear carbon chain, as in Fig. 2, calculated by various imaginary parts, under a finite bias voltage of 0.5 V, where the other calculation conditions are same as for Fig. 2b. The number pointed by the arrow denotes a grid spacing (eV) corresponding to the number of grid points. The reference values are obtained from calculations with a large number of grid points. (b) Mulliken populations in the carbon chain. The sequential numbers 1 and 12 correspond to the most left- and right-hand side atoms in the central region $C$, respectively.

Figure 4

(Color online) (a) Currents of the linear carbon chain calculated by the SCF calculations (solid line) and the interpolation scheme (dotted line). (b) Transmission of the linear carbon chain under a bias voltage of 0.3 V, calculated by the SCF calculations (solid line) and the interpolation scheme (dotted line). The imaginary part of 0.01 and the grid spacing of 0.01 eV are used for the integration of the nonequilibrium term in the density matrix. The other calculation conditions are the same as for Fig. 2b.

Figure 5

(Color online) 8-ZGNR with (a) a FM junction and (b) an AFM junction together with the spatial distribution of the spin density at the source-drain bias voltage $V_{\text{b}}=0\text{V}$. The zigzag edges are terminated by hydrogen atoms. The isosurface value of $\left|0.002\right|$ is used for drawing the spin density.

Figure 6

(Color online) $I-V_{\text{b}}$ curves for (a) the up-spin and (b) the down-spin states in the FM junction and (c) the up-spin and (d) the down-spin states in the AFM junction of 7- and 8-ZGNRs.

Figure 7

(Color online) Spin moments of the edge carbons in the central region $C$ for the FM junction of 8-ZGNR under various applied source-drain bias voltages $V_{\text{b}}$.

Figure 8

(Color online) (a) Calculation model for a tunneling junction consisting of four MgO(100) layers, being the $C_0$ region, sandwiched by Fe(100), being the $L_0$ and $R_0$ regions. The red large and small and blue circles denote Fe, O, and Mg atoms, respectively. (b) The net charge and spin magnetic moment of a Fe or Mg atom, belonging to each layer in the parallel magnetic configuration between the left and right leads. The position in the horizontal axis exactly corresponds to that of the layer in the above figure.

Figure 9

(Color online) $\mathbf{k}$-resolved transmission at the chemical potential for (a) the majority-spin state of the parallel configuration, (b) the minority-spin state of the parallel configuration, and (c) a spin state of the antiparallel configuration, respectively. For the calculations, $\mathbf{k}$ points of $120\times 120$ were used.

Figure 10

(Color online) (a) Optimized geometric structure of $\left({\text{LaMnO}}_3\right)/\left({\text{SrMnO}}_3\right)$ superlattice with the ferromagnetic spin configuration of Mn sites. The middle-sized blue and small red circles denote Mn and O atoms, while Sr and La atoms are denoted by the mark. DOS of the $\left({\text{LaMnO}}_3\right)/\left({\text{SrMnO}}_3\right)$ superlattice for (b) ferromagnetic and (c) A-type antiferromagnetic configurations.

Figure 11

(Color online) Speed-up ratio in the parallel computation of the calculation of the density matrix for the FM junction of 8-ZGNR by a hybrid scheme using MPI and OPENMP. The speed-up ratio is defined by $T_1/T_p$, where $T_1$ and $T_p$ are the elapsed times by a single core and a parallel calculations. The cores used in the MPI and OPENMP parallelizations are called process and thread, respectively. The parallel calculations were performed on a Cray XT5 machine consisting of AMD opteron quad core processors (2.3 GHz). In the benchmark calculations, the number of processes is taken to be equivalent to that of processors. Therefore, in the parallelization using 1 or 2 threads, 3 or 2 cores are idle in a quad core processor.