One-dimensional materials such as nanowires and nanotubes, which have unique electronic properties due to the quantum confinement effect, are expected to be applied to electronic and spintronic devices, optoelectronic circuits, and biosensors.Kong et al. (2000); Duan et al. (2001); Kind et al. (2002); Li et al. (2006); Patolsky et al. (2006); Liao et al. (2006); Ivanovskaya et al. (2007) In recent years, large-scale electron-transport calculations have been indispensable for designing the functionality of electronic devices. Although first-principles calculations based on the density functional theory (DFT)Kohn (1999); *PhysRev136-00B864; *PhysRev140-0A1133 allow us to accurately evaluate electron-transport properties of atomic systems, their target is typically limited to small systems because of the heavy computational cost arising from calculating the Green’s functions, which is intrinsic cubic scaling with system size. Limited-scale system is sufficient to simulate simple systems, such as locally perturbed bulk regions or interfaces. However, it is not suitable for simulating complex systems, such as bulk with a realistic defect density, amorphous structure, and interfaces with lattice mismatches between different crystalline materials. Thus far, to circumvent this restriction, transport calculations for large systems containing several thousands of atoms have been performed within atomic-basis formalism and tight-binding (TB) formalism based on DFT.Gunst et al. (2017); Ducry et al. (2019); Calogero et al. (2019); Jippo et al. (2016) Since the Hamiltonian matrix is dense in these formalisms, the computational cost of the inverse matrix for obtaining the Green’s function matrix by the direct method is proportional to the cube of the matrix dimension. Moreover, there are unavoidable problems such as the incompleteness of basis sets and inefficiency in parallel computing. However, real-space finite-difference (RSFD) formalism is also recognized as suitable for large-scale calculations requiring high computational accuracy because it has a high affinity to massively parallel architectures and can avoid problems arising from the basis sets.Fujimoto and Hirose (2003); Hirose et al. (2005); Ono et al. (2012); Iwase et al. (2015); Egami et al. (2015); Ono and Tsukamoto (2016); Tsukamoto et al. (2017); *PRB97-115450; *PRB98-195422; Egami et al. (2019)

Fujimoto and Hirose developed the overbridging boundary matching (OBM) method based on the RSFD formalismFujimoto and Hirose (2003) by exploiting the advantages for electron-transport calculations, where a whole system is divided along the $z$ direction (see Fig. 1) into three regions: the left electrode, the transition region and the right electrode. The electron-transport calculation refers only to subsets of the Green’s function matrix of the transition region. In the OBM method, the subsets are calculated efficiently using a shifted conjugate-gradient (SCG) solverIwase et al. (2015) since the Hamiltonian matrix in the RSFD formalism is sparse. In the SCG solver, the computational cost for calculating the subsets is ${\cal{O}}(N_{E}N_{xy}^{3}N_{z})$ when the algorithm proposed by Takayama et al.Takayama et al. (2006) is adopted, where $N_{xy}$ and $N_{z}$ are the numbers of grid points in the $xy$ plane and the $z$ direction, respectively, and $N_{E}$ denotes the number of energy values $E$ to be treated at once by the SCG method. For systems with a large $N_{z}$, the computational cost for matrix-vector operations, which is ${\cal O}(N_{xy}^{3}N_{z}^{2})$, becomes dominant.¹¹1Although the computational cost of the subsets of the Green’s function matrix is generally dominant, the cost of coefficient matrix-vector operations dominates in the case that $E$ becomes close to the eigenvalue of the truncated Hamiltonian matrix. The matrix-vector operation cost also becomes dominant when $N_{z}$ is large because the density of the eigenvalues in energy increases. Although the maximum order of the computational cost is reduced from ${\cal O}(N_{xy}^{3}N_{z}^{3})$ to ${\cal O}(N_{xy}^{3}N_{z}^{2})$ by the OBM method, the transport calculation has room for further reducing the computational cost to handle more realistic and longer one-dimensional systems.

In this letter, we propose an efficient computational procedure based on RSFD formalism to evaluate the electron-transport properties of long systems containing more than 100,000 atoms without accuracy deterioration. We verify the computational accuracy and performance of the proposed method. Additionally, to exemplify the efficiency of the proposed procedure we demonstrate large-scale electron-transport calculations for double-walled carbon nanotubes (DWCNTs) composed of 196,608 atoms, which make, to the best of our knowledge, the largest system in the first-principles electron-transport calculation.

Thouless and Kirkpatrick proposed an efficient calculation of subsets of a Green’s function matrix for a transition region within TB formalism as the recursive Green’s function (RGF) method.Thouless and Kirkpatrick (1981) The RGF method constructs a transition region by arranging multiple parts along the $z$ direction as illustrated in Fig. 1, and calculates only necessary subsets by recursively combining the Green’s function matrices of the adjoining small parts one by one. Since the calculation of the Green’s function matrices for the small parts is less costly than that for the extended transition region, the RGF method is capable of handling large systems and is widely used for calculating the transport properties of atomic-basis models and TB models.Lewenkopf and Mucciolo (2013); *JComputPhys215-000741; Ferry and Goodnick (1997); Sols (1995); Metalidis (2007); Souza et al. (2012); Shin et al. (2016); Rojas et al. (2019)

Now, we discuss the advantage of the RGF method in transport calculations within RSFD formalism. The transport properties are evaluated using the Green’s function of the transition region suspended between semi-infinite electrodes, $\mathbf{G}_{\mathrm{T}}$, which is constructed by the Green’s function associated with the truncated transition region, $\mathbf{\cal G}_{\mathrm{T}}$, and the self-energy terms of the electrodes.Ono et al. (2012) The Green’s function matrix $\mathbf{\cal G}_{\mathrm{T}}=(E\mathbf{I}_{\mathrm{T}}-\mathbf{H}_{\mathrm{T}})^{-1}$, in which $\mathbf{H}_{\mathrm{T}}$ is the truncated Hamiltonian matrix of the extended transition region consisting of $P$ parts (see Fig. 1), is expressed within the framework of RSFD formalism using the norm-conserving pseudopotentials (NCPPs) as

where $\mathbf{H}_{i}$ is the truncated Hamiltonian matrix for the transition region of the $i$th part. For simplicity, we assume that all $\mathbf{H}_{i}$ ($i=1,\ldots,P)$ have the same dimension, the interactions between the second- and further-neighboring parts are zero, and the non-zero matrix, $\mathbf{B}_{i}$, represents the interaction between the adjoining parts. We define the dimensions of $\mathbf{B}_{i}$ and $\mathbf{H}_{i}$ as $M_{i}\times N_{i}$ and $N_{xy}N_{z}^{(\text{part})}$, respectively, with $N_{z}^{(\text{part})}=N_{z}/P$; $\mathbf{I}_{\mathrm{T}}$ denotes an $N_{xy}N_{z}$-dimensional identity matrix.

Since the Hamiltonian matrix of RSFD formalism is highly sparse, it is more advantageous to calculate the subsets of Green’s function matrix of the transition region by iterative solvers than it is to compute the entire Green’s function by direct inversions of the Hamiltonian matrix. According to Ref. Ono et al., 2012, the electron-transport properties can be estimated using submatrices located at the four corners of the Green’s function matrix, $\mathbf{\cal G}_{\mathrm{T}}$. Now we will combine the two adjoining parts to derive the Green’s function matrix associated with a truncated Hamiltonian composed of the $i$th and ($i+1$)th parts, which is denoted by the $2N_{xy}N_{z}^{(\text{part})}$-dimensional matrix $\mathbf{\cal G}_{i:i+1}$. As described in Supplemental Material, the submatrices located at the upper-left, upper-right, lower-left, and lower-right corners of $\mathbf{\cal G}_{i:i+1}$ are given as

respectively. $\mathbf{I}_{M_{i}(N_{i})}$ denotes an $M_{i}(N_{i})$-dimensional identity matrix, and $\mathbf{\cal G}_{i:i}$ represents the $N_{xy}N_{z}^{(\text{part})}$-dimensional Green’s function matrix associated with $\mathbf{H}_{i}$, that is, $\mathbf{\cal G}_{i:i}=(E\mathbf{I}_{i}-\mathbf{H}_{i})^{-1}$, with $\mathbf{I}_{i}$ being an $N_{xy}N_{z}^{(\text{part})}$-dimensional identity matrix. As shown in Eqs. (2)–(5), we only need the four corners of $\mathbf{\cal G}_{i:i}$, which are efficiently obtained by using the SCG method. In most cases, $M_{i}=N_{i}=N_{\text{f}}N_{xy}$, with $N_{\text{f}}$ being the order of the finite-difference approximation.Chelikowsky et al. (1994a); *PRL72-001240 Repeating this procedure to connect the parts numbered from $i$ though to $j$, we can construct the four submatrices of a $(j-i+1)N_{xy}N_{z}^{(\text{part})}$-dimensional matrix $\mathbf{\cal G}_{i:j}$.

Now, we estimate the computational cost of the proposed method. We assume that $\mathbf{B}_{\mathrm{L(R)}}$ (see Fig. 1) and $\mathbf{B}_{i}$ for $i=1,\dots,P-1$ are $M$-dimensional square matrices with $M=N_{\text{f}}N_{xy}$. In the case of the RGF method based on the atomic-basis formalism, $\mathbf{H}_{\mathrm{T}}$ is dense; that is, $M$ is not sufficiently smaller than $N_{xy}N_{z}^{(\text{part})}$. Therefore, the iterative solvers are not more efficient than direct inversion, and there is not much benefit in using Eqs. (2)–(5) to reduce the computational cost. Furthermore, the overlap matrix is generally not any scalar matrix; hence, the SCG method cannot be used to solve the linear systems. However, the RGF algorithm must be compatible with RSFD formalism since $\mathbf{H}_{i}$ is a block-diagonal matrix and a sparse matrix.Hirose et al. (2005); Chelikowsky et al. (1994a, b) The proposed method can be implemented without negatively affecting the benefit of the SCG method, in which the matrix-vector operations for shifted energy points are omitted, and only the matrix elements of subsets are updated for shifted energy points during SCG iteration. The computational cost of the matrix-vector products with a sparse matrix increases linearly with $N_{z}$, and the number of SCG iterations is proportional to $N_{z}$. Therefore, the computational cost for the extended transition region is reduced from ${\cal O}(N_{xy}^{3}N_{z}^{2}\bigr{)}={\cal O}(N_{xy}^{3}(N_{z}^{\text{(part)}}P)^{2}\bigr{)}$ to ${\cal O}(N_{xy}^{3}(N_{z}^{\text{(part)}})^{2}P\bigr{)}$ using the RGF and the SCG methods since the cost is ${\cal O}(N_{xy}^{3}(N_{z}^{\text{(part)}})^{2}\bigr{)}$ for each of the $P$ parts; that is, the increase in computational cost can be suppressed linearly for $P$. Although the computational cost for calculating Eqs. (2)–(5) is ${\cal O}(PM^{3})$, where $M=N_{\text{f}}N_{xy}$, with $N_{\text{f}}$ being sufficiently smaller than $N_{z}^{(\text{part})}$. Consequently, further computational cost reduction for $N_{z}$ is achieved by making use of the advantages of the RGF method and SCG method within the RSFD formalism in comparison with the conventional procedure.Ono et al. (2012)

We evaluated the computational error between the conductance spectra using the Green’s function matrices obtained by the proposed method and the conventional methodOno et al. (2012) for the BN-doped (4,4)@(8,8) DWCNT to verify that the proposed method does not deteriorate the computational accuracy. The unit cell contains 192 atoms, and boron and nitrogen atoms are substitutionally co-doped into the outer tube along the circumferential direction. The dimensions are $L_{x(y)}=19.5~{}{\text{\AA}}$ and $L_{z}^{\text{(part)}}=9.84~{}{\text{\AA}}$. For this verification, a supercell has twice the dimension in the $z$ direction that of the unit cell, that is, $L_{z}^{\text{(ext)}}=2L_{z}^{\text{(part)}}$. Using the electronic structure calculation code rspaceHirose et al. (2005); Ono and Hirose (2005) based on RSFD formalism within the DFT framework, the effective potential of the unit cell is calculated with $1\times 1\times 10$ $k$-points, $N_{xy}=80\times 80$, $N_{z}^{(\text{part})}=40$, $N_{\text{f}}=4$, the NCPPs,Kobayashi (1999); Troullier and Martins (1991); *PRL48-001425 and the local-density approximation (LDA).Vosko et al. (1980) The Green’s function calculation is collectively performed using the SCG method with $N_{E}=41$. Figure 2 shows the conductance spectra of the supercell evaluated using the proposed method with $P=2$ and the conventional method. There are no notable errors more than $2.5\times 10^{-7}~{}{\text{G}_{0}}$.

Finally, we compared the computational times in calculating the Green’s function submatrices for several systems using the proposed method with those using the conventional method. Here the computational times are measured for systems in which (001) Si bulk is doped with B atoms. We consider a primitive unit cell of (001) Si bulk, which contains four Si atoms and has the dimensions of $L_{x(y)}^{\text{(part)}}=3.84~{}{\text{\AA}}$ and $L_{z}^{\text{(part)}}=5.43~{}{\text{\AA}}$. We measured the CPU time for calculating the submatrices of the Green’s function matrix in a supercell with $L_{z}^{\text{(ext)}}=PL_{z}^{\text{(part)}}$, in which the $P$ unit cells are arranged along the $z$ direction, ($P=1,2,4$ and 8), and any of the four Si atoms in each unit cell is randomly replaced with a B atom. The effective potential of the unit cell is calculated using the rspace code with $8\times 8\times 12$ $k$-points, $N_{xy}=12\times 12$, $N_{z}^{(\text{part})}=16$, $N_{\text{f}}=4$, the NCPPs, and the LDA. We prepared small parts with $L_{x(y)}=n_{x(y)}\times L_{x(y)}^{\text{(part)}}$ by duplicating the unit cell $n_{x(y)}$ times in the $x(y)$ direction with $n_{x}=n_{y}=1,2$, and 4 to verify the influence of the cross-sectional size on the CPU time. As summarized in Table 1, the CPU time for calculating the four submatrices of $\mathcal{G}_{1:P}$ with $N_{E}=16$ using the proposed method increases in proportion to $P$. We confirmed that the proposed method enabled us to notably reduce the computational cost and avoid computational difficulties in handling long systems. Furthermore, one can see that the time ratio is not attenuated even in large cross-sectional systems. Since the proposed method has high weak-scaling efficiency for increasing $P$, the benefit of this method becomes more remarkable as the system size increases.

With the above verifications, it is confirmed that electron-transport properties can be estimated efficiently without the deterioration of computational accuracy using the proposed method. As a practical demonstration of the method, we investigated the electron-transport properties of DWCNTs depending on the tube length. As illustrated in Fig. 3, eight DWCNT unit-cell models with different relative positions of the two BN dimers are prepared. Here, the transport properties are evaluated for a model in which the same unit cells (Model 4) are periodically combined (periodic model) and in which different unit cells (Model 1–8) are randomly arranged (random model). For the periodic model, the transport calculations are performed for the DWCNTs with $P=1$–$1,024$ to investigate the changes in the conductance spectrum according to the number of unit cells $P$ in the transition region. The numerical conditions and procedures to obtain the effective potential of each unit cell are the same as those mentioned above. Here, we ignore the electron-phonon couplings and evaluate the transport properties described in the ballistic regime.Ishii et al. (2010)

As shown in Fig. 4, for systems with small $P$, the conductance quantization is observed at the low incident energy, where there are four conduction channels with a transmission probability of almost unity. In the low energy region, the incident electrons flow through the intra-states of the inner and outer tubes, as depicted in Fig. 5(a). Therefore, there are high-transmission channels due to the inner-tube states that are insensitive to the scattering potential created by the BN dimers in the outer tube. However, the channel transmissions are notably reduced at the high incident energy because of the inter-states between the inner and outer tubes, as illustrated in Fig. 5(b). As $P$ increases, dips appear in the conductance spectra even at the low incident energy, and they gradually become prominent. This reflects the energy gap in the energy band structure that appears in an infinite DWCNT with a periodic scattering potential, as in the Kronig-Penney model. On the contrary, for the random model with $P=1,024$ (196,608 atoms, $\sim{1.0}~{}{\text{$\mu${m}}}$), overall conductance is suppressed as plotted in Fig. 4. Although the transmission in most conduction channels is notably reduced due to the random scattering potential component widely distributed in DWCNT, there are still channels with high transmission to which the intra-states of the inner tube contribute. Here, the random model with $P=1,024$ is prepared by randomly combining four different random models with $P=64$ and 96. We also performed calculations using random models with different configurations, but the results are essentially the same. Performing the large-scale electron-transport calculation, we can discover this difference in the conductance spectra between the randomly and periodically doped models.

In summary, we developed an efficient first-principles algorithm for evaluating the electron-transport properties of long systems consisting of a vast number of atoms. The Green’s function of the whole transition region extended towards the transport direction was obtained by recursively combining the Green’s functions of the adjoining parts one by one. The computational cost for calculating the submatrices of the Green’s functions required to estimate the transport properties is proportional to the number of the combined parts. As a practical application, we demonstrated the electron-transport calculations of BN-doped DWCNTs containing up to 196,608 atoms. The proposed method develops our understanding of experimental studies on the research and development of devices using low-dimensional materials, as reported by Refs. Franklin and Chen, 2010; Franklin et al., 2010; Qiu et al., 2017; Sun et al., 2019; Schmidt et al., 2019. Moreover, even in a submicron scale structure containing millions of atoms, it is possible to efficiently evaluate the electron-transport properties without any deterioration of computational accuracy.