Simple and efficient way of speeding up transmission calculations with k-point sampling

Introduction

Calculations of electronic conductance based on first principle methods such as density functional theory (DFT) provide a valuable tool in order to gain insights into electronic transport in nano-conductors and comparison to experiments without employing fitting parameters. This is for example the case in the field of single-molecular devices [1]. Popular methods are based on DFT in combination with the non-equilibrium Green’s function approach (DFT-NEGF), see, e.g., [2-4], or scattering wave-function approaches [5]. The electrodes in such calculations are typically treated employing periodic boundary conditions in the direction transverse to the transport direction with a corresponding k-point average of the electronic states and transmissions. This means that for each transverse k-point the system essentially behaves as a one-dimensional conductor with diverging density of states and discontinuities in the transmission function at energies corresponding to band on-sets/channel openings. It is well known that often in order to obtain smooth, well-converged density of states and transmissions as a function of energy, a substantial number of transverse k-points are needed due to the rapid variations of these functions for individual k-points [6]. Certain quantities, for example the Seebeck coefficient and thermo-electric figure of merit (ZT), are based on the detailed behavior of the transmission [7,8] and thus exceedingly sensitive to energy and k-resolution of the calculations. This can amount to a significant computational burden for large systems treated by first principle methods. Thus it is highly interesting to devise simple ways to make this more efficient and get maximum information from the data.

Sophisticated methods to tackle this include the transformation to a smaller basis-set using maximally localized Wannier functions [9], or to construct a optimized minimal basis-set and using this to determine the transmission [10]. Both require one to examine the details of the chemical bonds in the system, relevant energy windows, and storing wavefunctions, which tend to be elaborate. In this paper we present a simple and efficient post-processing interpolation scheme which can significantly speed up the convergence with respect to k-points. We illustrate the method by applying it to various graphene-based nano-structures which are prone to bad convergence due to its vanishing density of states at the Fermi level.

In the remaining parts of the paper we first explain the workings of the interpolation scheme in section Results and Discussion, while we investigate various test cases in section Example cases, and finally discuss limitations to the scheme and conclude in section Conclusion.

Results and Discussion

Description of the method

The use of computationally “expensive” first principles DFT-NEGF calculations for determining the transmission through nano-structured systems is limited by the amount of time one can afford to spend on the k-grid resolution. Often the rough behavior of the transmission is already seen with a limited number of k-points but the convergence of the average is slow since the functions are changing abruptly with energy, e.g., around a band onset or a resonance. The position of the abrupt feature will typically shift in a smooth way with changing k-point, but a linear interpolation of the curves between two consecutive k-points will be of little use since it will simply contain, say, half of each abrupt feature. Instead, we propose an interpolation scheme which can make use of a coarse, non-converged k-grid, and thus reduce the computation cost simply by using a “clever” technique to approximate the transmission curves for intermediate k-points. The method does not magically guess the correct interpolated curves, and one has to have a reasonable amount of k-point resolved transmission curves in order to obtain a useful result, but smooth averaged curves can be obtained, as we will illustrate below, using a significantly smaller number of k-points. The interpolation is done by using a shortest-path solver to determine a correspondence between two k-adjacent curves. The correspondence is then used to find intermediate curves which can be used to determine the averaged transmission. The proposed interpolation scheme consists of three separate steps, which will be described in detail in the following sections.

In order to show the validity of our proposed scheme we have determined the transmission through pristine graphene for increasing number of transverse k-points (see Figure 1a). The computational details are given in section Example cases. We compare N = 6, 8, …, 54 to a well-converged calculation with N = 600 k-points. By applying our algorithm we significantly reduce the mean absolute deviation from a fully converged transmission, as shown in Figure 1b. We see that the interpolated transmissions converge much quicker than the raw data. Fitting data with a power law, we can estimate the amount of k-points needed in order to obtain a deviation of less than one percent, thus giving a speed-up factor of η = 220/41 = 5.37.

Transform of data

In the following we will outline the method in general terms and denote the data points by (x, y), corresponding to the transmission function data point, (E, T_k(E)), for given k-point and energy in the concrete examples. Initially, we transform the set of data points (x, y) into points with a maximum Euclidian distance . The new data points span line segments (when connected) on the curve consisting of N different individual points. In practice the transformation of data is done by looping through all segment lengths L and inserting additional data points if L>. Figure 2a shows data where the distance between points a, c is much larger than , and thus requires extra points between points a, c. The optimal value of the segment length depends on the given data. A value close to the median of the original point distance is a good starting guess, and is chosen as the default value.

The two curves in Figure 2a are required to smoothly interpolate into each other – a process which can be estimated by hand (gray arrows). The proper path is found by overall minimization of the interpolation distance, which will described in the following subsection.

Correspondence between curves

Given two data sets, i.e., the curves , , we construct a matrix containing the weighted Euclidean distances between points on opposite curves,

where i = 1 … N₁, j = 1 … N₂ denote all points on the curves , , respectively. Thus, the dimension of the matrix is N₁×N₂. The weights w_x,y are added to ensure a degree of tunability when interpolating. This is due to the fact that x and y are different quantities and thus have different scales. In order to interpolate one curve into the other we need to determine the shortest path for all points on either curve. This is done by considering the distance matrix as a landscape and moving from the start (1, 1) to the destination (N₁, N₂) using the smallest cost possible. This is essentially identical to finding the minimum energy path on an potential energy surface, which can be done using either string methods or the nudged elastic band approach [11]. We will instead use Dijkstra’s algorithm, which loops through coordinates in the distance landscape and iteratively compares all found paths until the target point is found. Each iteration is done while keeping track of the cost and path. For more information see [12]. The coordinates along the returned shortest path is saved in a correspondence matrix containing two columns (, ), which is used in the final step of our routine. The distance matrix in Figure 2b shows the shortest path as a white line propagating along the distance landscape minimum (black/dark blue).

Interpolation of curves

Once the data is transformed and a correspondence matrix is obtained we interpolate the curves , . Instead of a simple linear interpolation we use the correspondence matrix to obtain intermediate curves. A simple linear, discretized interpolation is used, i.e.,

where α is the interpolation parameter between zero (curve ) and one (curve ), and x_i[M] means that we use the indices M of the vector x_i. The found intermediate data sets (x_12α, y_12α) can then be sampled back to the original grid. In practice we loop over α to generate the different interpolated curves. Figure 2c shows the data in Figure 2a interpolated into N = 10 intermediate curves and transformed back onto the original grid.

We note that the choice of weights w_x, w_y depends on the input data scales. Changing the values can highly affect the outcome of the shortest-path solver, since the distance landscape is changed. Usually, it is advisable to rescale the data (using the weights) so that the two data ranges are comparable. Similarly, the length has to be chosen wisely: A large value can result in crude interpolations while a too small value makes the algorithm too time-consuming.

Finally, we note that the algorithm described in the previous subsections allows us to interpolate data in general. We can apply the algorithm specifically to transmissions and DOS by providing it with the needed data. In the case of TranSiesta and the utility TBTrans we need to weight the interpolated curves using k-grid weights and sum to obtain average transmission curves. In the case of 3D transport we have a 2D k-grid, which can be investigated using bilinear interpolation. In the following section we apply the interpolation scheme to three example cases.

Example cases

The usefulness of the presented algorithm is showcased by considering transmission calculations through the simulated nano-systems shown in Figure 3: (a) a pristine graphene sheet, (b) a graphene nano-constriction [13], and (c) hydrogenated kinked graphene [14]. The shown structures have minimal unit-cells in the transverse transmission direction due to periodic boundary conditions, and carbon atoms are shown in black while hydrogen atoms are shown in white.

The transmission through the graphene sheet in Figure 3a is calculated with the TranSiesta simulation suite, which utilizes a localized basis set. For the present work a DZP basis set is used in conjunction with a mesh cut-off of 300 Ry, a force tolerance of 0.02 eV/Å, and a Monkhorst–Pack grid of 24 × 5 ensuring absolute convergence. Exchange and correlation is described using the PBE GGA functional [15]. A minimal transverse unitcell is used due to periodic boundary conditions, while an energy window of ±10 eV around the Fermi energy is considered both for a crude transverse transmission k-grid (N_k = 20) and a fully converged k-grid (N_k = 600). The coarse transmission spectrum is interpolated using the presented interpolation scheme, as can be seen in the upper row in Figure 4. The upper window of each column describes the full transmission versus k-point. A good agreement is seen between the interpolated and the converged data sets. A few places the interpolation guesses incorrectly (for instance around E = −4 eV for k = 0.25 ). These occur since the input data is too coarse. However, despite the few differences between interpolated data and converged data, the averaged curves have a remarkable resemblance and the mean deviation is below 2% (Figure 1b).

The transmission through the graphene nano-constriction and kinked graphene shown in Figure 3b and Figure 3c have been extracted from the original datasets (see [13,14]), and are shown in the middle and lower windows in Figure 4. As with the pristine graphene example we see a remarkable resemblance to the converged data set. The discrepancies are negligible, which stems from the fact that a finer k-point sampling has been performed in the original data. In the three cases we obtain speed-ups of approximately 5, 6, and 8, for the pristine graphene sheet, the graphene nano-constriction, and the kinked graphene sheet, respectively. Thus, we have demonstrated that by applying a simple post-processing interpolation scheme we can speed up convergence of roughly an order of magnitude.

Conclusion

We have presented a simple post-processing interpolation scheme which speeds up transmission calculations on nano-structured materials. The algorithm uses a shortest-path solver to determine the optimal interpolation of a set of k-dependent transmission curves, which ultimately can be summed to obtain a smoothed average transmission.

We note that as a post-processing tool the algorithm relies on the quality of the original data. Since this data is used as a base for the interpolation any fluctuations will be present in the interpolated data and thus propagate to the final result. The present implementation of the shortest-path solver is based on Dijkstra’s algorithm which is stable but very slow [12]. A far quicker implementation would be to use a heuristic to guide the shortest-path search in the distance landscape, thus changing the solver to an industry standard algorithm known as an A^*-search [16]. However, due to the complexity of the input data the construction of such a heuristic is not a straight-forward task.

By considering three sample cases we have demonstrated that it can speed up calculations by roughly an order of magnitude. We have illustrated the method using electron transport through graphene nano-structures and k-point averaging of transmission functions. However, the method is generally applicable also to phonon transport and to other functions such as density of states or other types of interpolation parameters such as electrostatic gating etc.

Our interpolation scheme can easily be implemented in already existing code. We provide a sample MatLAB code (Supporting Information File 1) that can read and interpolate data obtained from TranSiesta and TBTrans [2], which are built on-top of the ab initio software package Siesta [17].