A large-scale first-principles quantum transport simulation method using plane waves

TL;DR

This paper introduces a large-scale quantum transport simulation method using plane waves, enabling the study of systems with thousands of atoms and providing insights into nanoscale device properties.

Contribution

The paper develops a novel plane wave-based quantum transport simulation method capable of handling thousands of atoms, surpassing previous atomic basis set limitations.

Findings

Simulated nanowires with about 4000 copper atoms.

Analyzed effects of shape and point defects on transport properties.

Demonstrated efficiency of parallel algorithms for large systems.

Abstract

As the characteristic lengths of advanced electronic devices are approaching the atomic scale, ab initio simulation method, with fully consideration of quantum mechanical effects, becomes essential to study the quantum transport phenomenon in them. However, current widely used non-equilibrium Green's function (NEGF) approach is based on atomic basis set, which usually can only study small system with less than 1000 atoms in practice. Here we present a large-scale quantum transport simulation method using plane waves basis, based on the previously developed plane wave approach (Phys. Rev. B 72, 045417). By applying several high-efficiency parallel algorithms, such as linear-scale ground-state density function theory (DFT) algorithm, folded spectrum method, and filtering technique, we demonstrate that our new method can simulate the system with several thousands of atoms. We also use this method to study several nanowires with about 4000 copper atoms, and show how the shape and point defect affect the transport properties of them. Such quantum simulation method will be useful to investigate and design nanoscale devices, especially the on-die interconnects.

Tables

Table 1. Table 1: Number of atoms ( N Cu subscript 𝑁 Cu N_{\text{Cu}} ), the radius ( r 𝑟 r ) and length ( l 𝑙 l ) of the “neck” region, and the average transmission coefficients at the Fermi level of the electrodes ( T 𝑇 T ) of the system I-IV.

system	$N_{\text{Cu}}$	$r$ ($\mathrm{nm}$)	$l$ ($\mathrm{nm}$)	$T$	$T/r^2$ (${\mathrm{nm}}^{-2}$)
I	3996	0.5	2.0	0.179	0.716
II	3660	0.5	5.0	0.182	0.728
III	4140	0.75	5.0	0.350	0.622
IV	4139	0.75	5.0	0.352	0.626

Equations

$$I=\frac{2e}{h}\int [f(E-{\mu }_{\text{R}})-f(E-{\mu }_{\text{L}})]\sum _n{T}_n(E)\text{differential}E$$

$${\psi }_{\text{sc}}(r)=\{\begin{array}{ll}{\varphi }_n^{\text{R}\ast }(r)+{\sum }_{m\ne n}{B}_m^{\text{R}}{\varphi }_m^{\text{R}}(r),& x\to \infty \\ {\sum }_m{A}_m^{\text{L}}{\varphi }_m^{\text{L}\ast }(r),& x\to -\infty \end{array}$$

$$T_n(E)=\frac{{\sum }_m|A_m^{\text{L}}{|}^2[\text{d}{E}_m^{\text{L}}(k)/\text{d}k]{|}_{k_m^{\text{L}}}}{[\text{d}{E}_n^{\text{R}}(k)/\text{d}k]{|}_{k_n^{\text{R}}}}$$

$$R_n(E)=\frac{{\sum }_{m\ne n}|B_m^{\text{R}}{|}^2[\text{d}{E}_m^{\text{R}}(k)/\text{d}k]{|}_{k_m^{\text{R}}}}{[\text{d}{E}_n^{\text{R}}(k)/\text{d}k]{|}_{k_n^{\text{R}}}}$$

$$(\hat{H}-E){\psi }_l(r)=w_l(r)$$

$$(\hat{H}-E{)}^2{\psi }_l(r)=(\hat{H}-E)w_l(r)=w_l^{\prime }(r)$$

$$F=\frac{1}{2}\text{expectationvalue}(\hat{H}-E{)}^2{\psi }_l-\text{innerproduct}{\psi }_l{w}_l$$

$${\psi }_l=\sum _i\frac{\text{innerproduct}{\varphi }_i{w}_l^{\prime }}{({\epsilon }_i-E{)}^2}{\varphi }_i$$

$$\begin{array}{rl}{M}^{\prime }& =M+\sum _{i=1}^N\left[\frac{1}{({\epsilon }_i-E{)}^2+{\Delta }_i^2}-M\right]\text{outerproduct}{\varphi }_i{\varphi }_i\\ & =M\left(1-\sum _{i=1}^N\text{outerproduct}{\varphi }_i{\varphi }_i\right)+\sum _{i=1}^N\frac{\text{outerproduct}{\varphi }_i{\varphi }_i}{({\epsilon }_i-E{)}^2+{\Delta }_i^2}\end{array}$$

$$\hat{H}(k_{1,2})|{\psi }_i(k_{1,2})⟩={\epsilon }_i(k_{1,2})|{\psi }_i(k_{1,2})⟩$$

$$\hat{H}(k)\approx \alpha \hat{H}(k_1)+(1-\alpha )\hat{H}(k_2)$$

$$H_{ij}(k)=\text{matrixelement}{\psi }_i(k_1)\hat{H}(k){\psi }_j(k_1)$$

$$=\alpha {\epsilon }_i(k_1){\delta }_{ij}+(1-\alpha )\text{matrixelement}{\psi }_i(k_1)\hat{H}(k_2){\psi }_j(k_1)$$

$$H_{ij}(k)=\alpha {\epsilon }_i(k_1){\delta }_{ij}+(1-\alpha )\sum _m{U}_{im}^{\ast }{U}_{jm}{\epsilon }_m(k_2)$$

$${\psi }_i(k)=\alpha \sum _j{c}_j{\psi }_j(k_1)+(1-\alpha )\sum _{j,k}{c}_j{U}_{jk}{\psi }_k(k_2)$$

$${\rho }_E(r)=\sum _{i=1}^N|{\psi }_{\text{sc}}(r){|}^2$$

Full text

quantum transport, ab initio simulation, plane waves basis, linear-scale

I INTRODUCTION

The characteristic scale of advanced electronic devices have reached nanoscale. As the International Roadmap for Devices and Systems (IRDS)IRD shows, the active channel length of an electronic device has been downscale to tens of nanometers. The traditional technology computer aided design (TCAD) method is based on continuum level models, which cannot be applied to such small systems where the atomistic quantum effects are dominate. Instead, atomistic modeling and simulation method based on quantum mechanics become necessary.

In the transport simulation, a two-probe model, which contains two semi-infinity electrodes and a central region for scattering, is widely used. The current in this model can be calculated by using the Landauer-Büttiker formalism

[TABLE]

Here $f$ is the Fermi–Dirac distribution function, $\mu_{\text{L}}$ and $\mu_{\text{R}}$ are the Fermi level of left and right electrodes, and $T_{n}(E)$ is the transmission coefficient for the $n$-th right electrode band at energy $E$. Since $\mu_{\text{L/R}}$ can be obtained by using ab initio bulk calculation, the key problem of transport simulation is to calculate $T_{n}(E)$. Currently, method based on non-equilibrium Green’s function (NEGF) theoryBrandbyge et al. (2002) is the most widely used approach. Tight-binding model based NEGF approach can simulate a system containing more than 100000 atoms,Klimeck et al. (2007a, b) but it ignores the charge density self-consistent effect which is important at atomistic level. Furthermore, NEGF approach combining with ab initio density functional theory (DFT) prefer a localized atomic basis set rather than plane waves. However, a localized basis is less accurate and less flexible than a plan wave one in DFT calculations.Davidson and Feller (1986) Although there are several approaches to generate localized Wannier functions from the plane wave calculation, many approximations must be made and many of these method are computationally expensive.Calzolari et al. (2004) The widely used DFT-NEGF approach usually can only simulate system less than 1000 atoms, which is not enough for a real nanoscale device.

Given these circumstances, it will be interesting to test different methods. One approach is to directly solve the scattering states using plane wave basis and auxiliary periodic boundary conditions.Wang (2005); Garcia-Lekue and Wang (2006); Garcia-Lekue et al. (2015) This approach can be easily applied to the existing DFT code, and the computational resource it takes is similar to that of a standard ground-state calculation.Wang (2005) However, since the computational complexity of common DFT method is $O(N^{3})$ (here $N$ is the scale of system), and many eigenstates of the system need to be calculated, it can only be used to study the system contains hundreds of atoms in practice. In this work, we combine this plane waves simulation method with linear-scale ground-state DFT algorithm, as well as the folded spectrum methodWang and Zunger (1994) and filtering techniques, to solve problems with thousands of atoms. In particular, we study on-die interconnect of copper nanowires, with several thousands atoms, and show the shape and point defect effect in them. Our new method will be useful to study large-scale nano-electronic devices, especially for nanoscale interconnects.

II COMPUTATIONAL PROCEDURE

The system we study in this work is copper nanowires, which are commonly used for on-chip interconnects of nowadays integrated circuit (IC). And such interconnect is the main source of heat generation, is also facing difficulty for further downscaling. In this work we use four systems with the same electrodes(which is a copper nanowire with radius of $1\text{\,}\mathrm{n}\mathrm{m}$, and its cross sections are shown in FIG. 1(a), (b)) to study. The structure of four systems (which are called system I, II, III, and IV) are shown in FIG. 1(c)-(f), the transport direction are all along the $[100]$ direction of bulk copper, which is the $x$-direction in these structures. The main difference between system I, II, and III is radius and length of narrow part (we call it “neck”) in the central region, while system IV is created by introducing a vacancy defect in system III (as the arrow in FIG. 1(f) indicate). The number of atoms in these four systems are 3996, 3660, 4140, and 4139, respectively. The good to use these system is to investigate how the central part shapes and point defect affect the interconnect transmission.

In all the calcultions of this work, the local density approximation (LDA) with the Perdew-Zunger parameterizationPerdew and Zunger (1981) of the exchange-correlation functional and Fritz-Haber-Institute (FHI) norm-conserving pseudopotentialsFuchs and Scheffler (1999) are used, and the cutoff energy is set to $40\text{\,}\mathrm{R}\mathrm{y}\mathrm{d}$ ($\sim$540\text{,}\mathrm{e}\mathrm{V}$$). All the calculations are done on the Titan supercomputerTia , while the central processing unit (CPU) and graphics processing unit (GPU) are AMD Opteron 6274 and Nvidia K20X, respectively.

II.1 The original approach

The central idea of the plane wave approach is to get scattering state wave funtions $\psi_{\text{sc}}$ by solving the Schrödinger equation $\hat{H}\psi_{\text{sc}}=E\psi_{\text{sc}}$ ($E$ is the given energy), with the following boundary conditions

[TABLE]

where $\phi_{n}^{\text{L/R}}({\bf\it r})$ and $\phi_{n}^{\text{L/R}*}({\bf\it r})$ are the left- and right-going running wave states in the left/right electrodes with energy $E$. The transmission and reflection coefficient can be calculated by

[TABLE]

To solve the scattering states, we introduce some “system states” $\psi_{l}$ which satisfy the equation

[TABLE]

where $w_{l}$ is a perturbation function that is nonzero only near the boundary far from the central region. The $\psi_{l}$ are decomposed to left and right electrode running wave states $\phi_{m}^{\text{L}}$ and $\phi_{m}^{\text{R}}$ at one left electrode unit cell $\Omega_{\text{L}}$ and one right electrode unit cell $\Omega_{\text{R}}$. Using such decomposition, we can then use linear combination of $\psi_{l}$ to generate the scattering states which satisfy the boundary conditions Eq. (2) (where the certain running wave coefficients are zero). We can then get $T_{n}(E)$ and thus solve the transport problem. More details of the original method can be found in Ref. Garcia-Lekue et al. (2015).

II.2 Electronic structure

The above approach can be easily implemented in the common DFT codes, and its computational time is similar to that of a standard ground-state calculation. A package called PEtot_transPWt has been developed using this approach based on PEtot codePEt . But common DFT algorithm has the computational complexity of $O(N^{3})$, that makes it very expensive for large systems. In practice, we can only calculate the electronic structure of a system with hundreds of atoms using the above approach. Fortunately, several linear-scale DFT algorithms have been developed.Kim et al. (1995); Goedecker and Colombo (1994); Li et al. (1993); Mauri et al. (1993); Ordejón et al. (1993); Yang (1991) Those algorithms have the computational complexity of $O(N)$, enable them for large system calculations. One such algorithm is called linear-scaling three-dimensional fragment (LS3DF) method, and has been used to study the electronic structure of the systems containing more than 10000 atoms based on plane wave expansions.Wang et al. (2008); Zhao et al. (2008); Wang (2014) This is benefited from the highly-scalable divide-and-conquer scheme implemented in the LS3DF code, which makes it not only linear scaling to the system size, but also linear scaling to the number of processors used in the computations.

In this work, we use a GPU version of LS3DF codeJia et al. (2017). Like most linear scaling DFT method, LS3DF is based on the nearsightedness of the DFT properties. Such nearsightedness is rigorous for semiconductors,Kohn (1996) but it is not clear whether LS3DF will still work for metallic system. We thus need first to check whether LS3DF works for the metallic system. Here we use a relatively small system of 488 copper atoms, whose structure is similar to those we study, to test the validity of LS3DF and its speed. We find that it takes about 500 minutes if we use the common Kohn-Sham DFT code (PEtot) with 360 CPU processors (5760 cores), but only 30 minutes if the GPU version of LS3DF code and 480 GPUs are used. In the LS3DF calculation, the system has been divided into 95 overlapping fragments, and each fragment contains $9\sim 334$ copper atoms without any surface passivation. To compare the results calculated by these two software, we calculate the average potential on the $yz$-plane ($\bar{V}(x)$) and show them in FIG. 2(a). We can find that the error between these two method are less than $0.02\text{\,}\mathrm{e}\mathrm{V}$ (as shown in FIG. 2(b)), which indicated that LS3DF gives almost the same result as PEtot and thus suggested that LS3DF works well for the metallic system.

In LS3DF calculation, a large system is first divided into small overlapping fragments. The wave functions and charge densities of individual fragments are calculated separately by specific groups of computer processor. These fragment charge densities are then patched together to yield the global charge density using a special patching scheme, The Poisson equation of the whole system is solved based on the patched global charge density. For a metallic system, a global Fermi energy is used to control the orbital occupations of all fragments. For a given accuracy, the size of the fragment will be fixed regardless of the total system size. The larger system just means larger number of fragments. Besides, larger number of computer processor means small number of fragments each computer processor group has to deal with. This leads to a double linear scaling property of LS3DF: linear both to the size of the system, and number of computer processor used. Due to the double linear scaling, LS3DF can calculate the large systems quickly. For instance, it only takes about 30 minutes to calculate the charge density of the central region of the system III (which contains about 3000 atoms) using 1920 GPUs, and the results are shown in FIG. 3.

II.3 The folded spectrum method and acceleration technique

After we solved the problem of LS3DF for self-consistent field solution of the charge density and potential of the system, we now turn our attention to solve the scattering state of the system. The preconditioned conjugate gradient (PCG) method can be quite efficient to solve the linear equation Eq. (5). Unfortunately, $(\hat{H}-E)$ is not positive definite, so the PCG method cannot be directly applied. In the original approach, we have calculated all the eigenstates of $\hat{H}$ below the given energy $E$ and depleted those eigenstates from $\hat{H}$ to make $(\hat{H}-E)$ positive definite. However, this approach only works for small system where all eigenstates below $E$ can be solved using iteration methods. For a system contains thousands of atoms, there could be tens of thousands of eigenstates below $E$ (which is usually near the Fermi level of electrodes), so it will be very expensive to calculate all of them. This becomes the main limitation to solve the scattering states using plane wave basis set.

To solve this problem, we use the folded spectrum methodWang and Zunger (1994), which apply another $(\hat{H}-E)$ on both sides of Eq. (5), then it becomes

[TABLE]

Since $(\hat{H}-E)^{2}$ is always positive definite, this equation can be directly solved by using the PCG method. For the PCG method, we use PCG to minimize the following quantity:

[TABLE]

For some system, such as a 2000 atoms black phosphorus nanosheet, the PCG iteration converges in several thousands steps when solving Eq. (6), thus it can be directly used. But in many cases, such as the copper system we study here, the PCG iteration is very hard to converge. This is due to the poor condition number of operator $(\hat{H}-E)^{2}$, especially when there are eigenstates with eigenenergy $\varepsilon_{i}$ close to $E$. For the four copper systems we study, more than 30000 steps of PCG is required to converge (which means the error, $\norm{(\hat{H}-E)\psi_{l}-w_{l}}$, is less than $1\text{\times}{10}^{-8}$, in atomic unit). This is rather expensive. Especially, since there is no way to massively parallelize the PCG iterations, the calculation takes a long time, without the benefits of massive supercomputer.

The solution of Eq. (6) can be expanded by all the eigenstates of $\hat{H}$

[TABLE]

Here $\phi_{i}$ and $\varepsilon_{i}$ are the $i$-th eigenstate and eigenvalue of $\hat{H}$. Thus the eigenstates around the given energy $E$ contribute the most to $\psi_{l}$. If we have these eigenstates, they can be used to accelerate the PCG convergency (for instance, they can be depleted from the Hamiltonian). In this case, those eigenstates are just used to accelerate the convergency, there is no need to calculate all the eigenstates below $E$, thus we need much smaller number of eigenstates.

To get a limited number of eigenstates around $E$, there are several methods one can use, one is the folded spectrum methodWang and Zunger (1994). However, we meet the same converge problem as solving Eq. (6). Here we adopted the Chebyshev filter diagonalization (ChebFD) algorithmWeiße et al. (2006); Pieper et al. (2016), to find the eigenstates within a energy window around $E$. In the ChebFD algorithm, a polynomial of the Hamiltonian $\hat{H}$, which is constructed by the linear combination of Chebyshev polynomials, is used to imitate a window function surrounding $E$, and this polynomial $p(\hat{H})$ can be viewed as a filter. Then this filter can be applied to a group of random wave functions $\{\phi_{i}^{(0)}\}$. Assume $\phi_{i}^{(0)}=\sum_{j}c_{ij}\phi_{j}$, where $\phi_{j}$ is the eigenstate of $\hat{H}$ with eigenenergy $\varepsilon_{j}$, then $\phi_{i}^{\prime(0)}=p(\hat{H})\phi_{i}^{(0)}=\sum_{j}c_{ij}p(\varepsilon_{j})\phi_{j}$, thus it will suppress the amplitudes of eigenstates with their $\varepsilon_{j}$ outside the window, while keeping the eigenstates inside the window. One can then use $\{\phi_{i}^{\prime(0)}\}$ to do a subspace diagonalization of $\matrixelement{\phi_{i}^{\prime(0)}}{\hat{H}}{\phi_{j}^{\prime(0)}}$, yield eigenstates $\{\phi_{i}^{(1)}\}$ of this subspace. This is one ChebFD step, and one can repeat this ChebFD step by applying $p(\hat{H})$ again to $\{\phi_{i}^{(1)}\}$. Several ChebFD iterations need to be used to get accurate eigenstates in practice. One major benefit of ChebFD algorithm is its parallelization. Because $p(\hat{H})\phi_{i}^{(n)}$ can be done independently for different $\phi_{i}^{(n)}$, the procedure is embarrassingly parallel. Although this algorithm may loss some eigenstates in the energy window, this leakage is not important for our case, since we only use them as preconditioner. FIG. 4 shows how the ChebFD algorithm works for system III, the largest one in our test systems. Here a 1500-degree of ChebFD filter (1500 order in $p(\hat{H})$ polynomial) and 1500 random initial wave functions are used. The energy window is set to $2.0\text{\,}\mathrm{e}\mathrm{V}$. As shown in FIG. 4, the overall error of the eigenstates decrease by about one order of magnitude with each step of ChebFD. Here we use 43200 CPU cores on Titan, and it takes about half an hour to run one step of ChebFD. In doing so, we have used 432 CPU cores for each $p(\hat{H})$. and the total number of processors have been divide into 100 groups. In this part of calculations, only CPU codes are available. We expect future implementation on GPU nodes can significantly accelerate this calculation, and make such calculation more routine.

When we get the eigenstates around $E$, we can use them to construct a projector and than apply it to Eq. (6) (deplete out such eigenstates), as we did in the original approach. But that require all the eigenstates are fully converged. Here, instead, we use them as preconditioner. As we know, a good preconditioner can significantly improve the convergency of PCG iterations. We have already used one $G$-space diagonal preconditioner. but it is only effective on the high energy plane wave at the higher end of $(\hat{H}-E)^{2}$. Now we have some eigenstates around $E$, we can use them to construct a much better preconditioner (which is closer to $(\hat{H}-E)^{-2}$)

[TABLE]

here $M(G)=\frac{1}{(G^{2}-E)^{2}+\Delta E^{2}}$ is our $G$-space diagonal preconditioner, and $M^{\prime}$ is our new preconditioner. Basically, in the subspace $\operatorname{span}\{\outerproduct{\phi_{1}}{\phi_{1}},\outerproduct{\phi_{2}}{\phi_{2}},\cdots,\outerproduct{\phi_{N}}{\phi_{N}}\}$, we use the $(\hat{H}-E)^{-2}$ inversion, while in the rest of the space, we use the $M$ preconditioner. This formalism looks almost the same as the projector one. But in the current approach, the eigenstates $\{\phi_{i}\}$ need not to be very accurate. In particular, we have $\Delta_{i}=\norm{\hat{H}\phi_{i}-\varepsilon_{i}\phi_{i}}$ is the error of the eigenstate $\phi_{i}$. In this way, we can use the approximated eigenstates to accelerate the PCG convergency.

Here we also use system III to test our accelerate technique, all the results are shown in FIG. 5. For that system, the linear equation error defined as $\norm{(\hat{H}-E)\psi_{l}-w_{l}}$ decreases very slowly if no eigenstates are used for acceleration, and it is still larger than $0.1$ (in atomic unit) after 10000 steps of PCG. Then we use eigenstates with different accuracy (which means generated by different steps of ChebFD, as shown in FIG. 4) to construct the preconditioner and then repeat the PCG iterations. Unfortunately, we find that eigenstates generated by 2 or 4 steps of ChebFD only accelerate the PCG convergency modestly, which suggests that eigenstates with low accuracy ($\Delta_{i}>$1\text{\times}{10}^{-5}) cannot be used for effective acceleration. Only when the eigenstates with $\Delta_{i}\leq$1\text{\times}{10}^{-6} are used in the calculations, and we find the iteration has been accelerated significantly. It takes about 6200, 4100, and 3500 PCG steps to suppress the linear equation error below $1\text{\times}{10}^{-8}$, when using eigenstates generated by 6, 8, and 10 steps of ChebFD, respectively. As shown in FIG. 4 , $\Delta_{i}$ of most eigenstates used here are less than $1\text{\times}{10}^{-6}$, $1\text{\times}{10}^{-8}$, and $1\text{\times}{10}^{-10}$ in these three cases. It takes only 2-3 seconds for one PCG step using 432 CPU cores, so the solving of the linear equation itself is not expensive if we already have the eigenstates. Our result also shows that the error of eigenstates ($1\text{\times}{10}^{-6}$ for this system) can be greater than the linear equation error ($1\text{\times}{10}^{-8}$ for this system), which will not be true in the projector depletion approach. This shows the advantage of using preconditioner approach.

II.4 Wave function interpolation technique

As discussed above, by using the LS3DF and folded spectrum method (with acceleration), we have solved two main problems to apply plane-wave transport method to large system. Nevertheless, there are still some detailed issues in practice worth to be discussed.

In a large system, the electronic structure of the electrodes can be quite complex. Since we need to calculate the running-wave states in the electrode, the band structure with many hundreds of $k$-points along the transport direction (which is $x$-direction in this work) is required. It turns out to be very expensive to do such calculations with so many $k$-points for an electrode with hundreds of atoms. Here we employ wave function interpolation to extract the eigenstates between two nearest $k$-points which have already been calculated. In this way, we can use relatively small number of $k$-points to get the accurate whole band structure. The details of the wave function interpolation technique are discussed in the following.

If we have already calculated the orthonormal eigenstates of two nearest $k$-points, $k_{1}$ and $k_{2}$, which means

[TABLE]

Since $k_{1}$ and $k_{2}$ are very close, so for a $k$ in $(k_{1},k_{2})$, the Hamiltonian $\hat{H}(k)$ can be expressed using a linear interpolation

[TABLE]

where $\alpha=\frac{k-k_{2}}{k_{1}-k_{2}}$. So the matrix elements of $\hat{H}(k)$ can be written as

[TABLE]

Here a transformation matrix $U$ between $\psi_{i}(k_{1})$ and $\psi_{j}(k_{2})$ can be used, which means $\psi_{i}{(k_{1})}=\sum_{j}U_{ij}\psi_{j}(k_{2})$. If all the eigenstates are included, $U$ will be the overlap matrix $U_{ij}=S_{ij}=\innerproduct{\psi_{j}(k_{2})}{\psi_{i}(k_{1})}$. In reality, only finite number of $\psi_{i}$ are used, thus the basis is incomplete, and $S$ is not a unitary matrix. A Gram-Schmidt process on $S$ is used to make a unitary $U$ out of $S$ matrix. Under this approximation, we can rewrite Eq. (12) as

[TABLE]

Then we can get the eigenstates $\psi_{i}(k)$ and $\varepsilon_{i}(k)$ by diagonalizing of $H(k)$, and $\psi_{i}(k)=\sum_{j}c_{j}\psi_{j}(k_{1})$. However, this formula will have a jump at $k_{2}$. To make it continuous, we can expand $\psi_{i}(k)$ by both $\psi_{i}(k_{1})$ and $\psi_{i}(k_{2})$, such as

[TABLE]

As we discuss above, all the system we study have the same structure of electrode, which contains 145 copper atoms per unit cell. In the ground-state calculation, we use 101 $k$-points along the $x$-direction and output all the eigenstates in real space. Then we use the above wave function interpolation technique, to add 9 $k$-points between every two nearest $k$-points. In this way, we get the band structure of 1001 $k$-points along the $x$-direction.

The band structure without and with wave function interpolation are shown in FIG. 6, which are represented by the blue dash and red solid line, respectively. It is easy to find that those two band structures are the same in most parts, but have some differences when it is blewd up. Overall, the bands with interpolation are much smoother than those without interpolation. We also find that there are some crossing bands in the original band structure, which become anitcrossing bands after interpolation. In the following, we will get running-wave states with higher accuracy using bands with interpolation.

Since the system we study are all metallic structures, the transport properties near the Fermi level of electrode are the most important. At the Fermi level (which has been adjusted to zero in FIG. 6), there are 52 running-wave states (their conjugate states are not counting), but this number changes to 50 after interpolation, as shown in FIG. 6(a). This is due to the detailed difference between those two band structures, especially the anticrossing bands, which is shown in FIG. 6(b). Each of these running wave states $\phi_{i}$ with energy of $E$ will be used to construct one $w_{i}$ by multiplying $\phi_{i}({\bf\it r})$ with a mask function $m(x)$, which is nonzero only near the artificial boundary of the system, far away from the central part of the interconnect.

III RESULTS AND DISCUSSION

Finally, we could use the calculated system states $\psi_{l}$ to construct scattering state $\psi_{\text{SC}}$ and thus solve the transport problem. Since all the three systems are metal, only their transport properties at the Fermi level of the electrodes are studied in this work. The transmission ($T_{n}$) and reflection coefficient ($R_{n}$) of each scattering states $\psi_{\text{sc}}(n)$ are shown in FIG. 7. Note we have made $T_{n}+R_{n}=1$.

As we discussed above, there are 50 left-going running wave states from the right electrode at the Fermi level (the black dots in FIG. 6(a)). FIG. 7 shows the transmission and reflection coefficient of the those incoming electrode states in the systems we studied. We find that there are always some states with relatively large $T_{n}$ (for example, the states with $k(x)=0.451$ and 0.475), while the magnitudes of $T_{n}$ of some states are greatly influenced by the “neck” structure of the system. Overall, system I and and II have similar transmission patterns than system III and IV, indicating the importance of the narrowest neck in controlling the transmission.

The left panel in FIG. 8 shows the scattering wave functions of the four interconnect from the same right-incoming wave (with $k(x)=0.583$) as indicated by the vertical arrow shown in FIG. 7. We find that the detailed structure in the central region (where the scattering occurs) influences the scattering wave functions significantly. System II, III, and IV have more components of wave function in the left electrode than that of system I, therefore their $T_{n}$ is much larger. We have also calculated the charge density of all the scatter states wave functions at the Fermi energy

[TABLE]

Here $N$ is the number of left-going running-wave states in the right electrode. These are shown in the right panel of FIG. 8. Since $\rho_{E}({\bf\it r})$ mainly distribute in the right side of all the systems we study, it suggests that the overall transmission coefficient of those system is not very large.

However, it is hard to find how the shape of central part influences $T_{n}$ for each individual scattering states, as we discussed above. But in general, system III and IV have larger $T_{n}$ than I and II. We can define the average transmission coefficient $T=\sum_{i=1}^{N}T_{i}/N$ to show this tendency. This $T$ also represents the overall conductivity of this interconnect. As shown in TABLE 1, the average transmission coefficient of system I and II are very close, which is about half of that of system III and IV. This $T$ increases with the “neck” diameter $r$. However, the $T/r^{2}$ ratio is larger for smaller neck interconnect, indicating the quantum mechanical effect, might enhance the transmission on top of the simple classical estimations. The other geometry parameter $l$, which measure the length of the “neck” region, is not so sensitive compare to the radius of the “neck” $r$ in terms of controlling $T$. $T$ only decreases slightly when $l$ increases from $2.0$ to $5.0\text{\,}\mathrm{n}\mathrm{m}$.

In system IV, we introduce a vacancy defect in system III, as shown in FIG. 8(d). We find that it only influences the scattering state of the same channel very slightly. Actually, it increases the overall $T$ slightly. It is rather surprising for such small influence given the small size of the neck. This phenomenon suggests that some point defect might not have detrimental effects on the transmission. All these phenomenon warrant more comprehensive future investigations.

IV CONCLUSION

In summary, we have improved a previously developed ab initio quantum transport simulation method based on plane wave basis set. By introducing the linear-scale ground-state DFT algorithm (LS3DF), folded spectrum method (with acceleration), and other techniques such as filtering and wave function interpolation, we are capable to simulate the transport properties of a system with several thousands of atoms. Since the computational complexity of our new method is nearly linear to the scale of system, the system scale will be able to be extended to tens of thousands with massively parallelize. The future implement on GPU nodes is also expected to accelerate the calculation.

By using this new method, we study the transport properties of several copper nanowire systems with about 4000 atoms, and show the shape and point defect effects in them. We find that the conductance of copper nanowire is mainly determined by the radius of narrow part, not the length, and some point vacancy defect might even improve the transmission. These simulations suggest the potential usage of our method for the study and design of nanoscale interconnect in the advanced electronic devices.

Ab initio simulation is useful to investigate and design the nano-electronic devices. Our method provide a new approach to study the transmission phenomenon in them. However, real semiconductor device, such as transistor, is more complex than then interconnect system we study here. Further development is expected to make our method can be applied to real devices with tens of thousands of atoms.

Acknowledgements.

L.W.W. was supported by the US Department of Energy, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division under contract no. DE-AC02-05-CH11231 within the beyond-Moore’s law LDRD project. M.Y., X.J., and S.S.L. were supported by China Key Research and Development Program (2018YFA0306101), National Natural Science Foundation of China (Grand Nos. 11774338, 11574304), Chinese Academy of Sciences-Peking University Pioneer Cooperation Team (CAS-PKU Pioneer Cooperation Team), the Youth Innovation Promotion Association CAS (Grand No. 2016109). This work used the resources of the National Energy Research Scientific Computing Center (NERSC) and Oak Ridge Leadership Computing Facility (OLCF) through the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) project NTI009.