GPUQT: An efficient linear-scaling quantum transport code fully implemented on graphics processing units

TL;DR

GPUQT is a GPU-accelerated quantum transport code that enables efficient, linear-scaling calculations of electronic properties in large disordered systems using advanced numerical techniques.

Contribution

The paper introduces GPUQT, a fully GPU-implemented quantum transport code that achieves linear scaling for large systems using novel numerical methods.

Findings

Achieves linear scaling in computation time and memory usage.

Successfully computes electronic transport properties in large disordered systems.

Demonstrates practical usage with example applications.

Abstract

We present GPUQT, a quantum transport code fully implemented on graphics processing units. Using this code, one can obtain intrinsic electronic transport properties of large systems described by a real-space tight-binding Hamiltonian together with one or more types of disorder. The DC Kubo conductivity is represented as a time integral of the velocity auto-correlation or a time derivative of the mean square displacement. Linear scaling (with respect to the total number of orbitals in the system) computation time and memory usage are achieved by using various numerical techniques, including sparse matrix-vector multiplication, random phase approximation of trace, Chebyshev expansion of quantum evolution operator, and kernel polynomial method for quantum resolution operator. We describe the inputs and outputs of GPUQT and give two examples to demonstrate its usage, paying attention to the interpretations of the results.

Equations

$${\sigma }^{\mathrm{KG}}(E)=\frac{2\pi \hslash e^2}{\Omega }\text{Tr}\left[V\delta (E-H)V\delta (E-H)\right],$$

$${\sigma }^{\mathrm{VAC}}(E,t)=e^2\rho (E){\int }_0^t{C}_{vv}(E,t)dt;$$

$$C_{vv}(E,t)=\frac{\frac{2}{\Omega }\text{Re}\left[\text{Tr}\left[U(t)V\delta (E-H)U(t{)}^{†}V\right]\right]}{\frac{2}{\Omega }\text{Tr}\left[\delta (E-H)\right]};$$

$$\rho (E)=\frac{2}{\Omega }\text{Tr}\left[\delta (E-H)\right],$$

$${\sigma }^{\mathrm{MSD}}(E,t)=e^2\rho (E)\frac{d}{2dt}\Delta X^2(E,t),$$

$$\Delta X^2(E,t)=\frac{\frac{2}{\Omega }\text{Tr}\left[[X,U(t){]}^{†}\delta (E-H)[X,U(t)]\right]}{\frac{2}{\Omega }\text{Tr}\left[\delta (E-H)\right]},$$

$$C_{vv}(E,t)=v_x^2(E)e^{-t/\tau (E)},$$

$${\sigma }_{sc}(E)=e^2\rho (E)v_x^2(E)\tau (E).$$

$$\lambda (E)=v_x(E)\tau (E).$$

$$D(E)=\frac{1}{2}\underset{t\to \infty }{\lim }\frac{d}{dt}\Delta X^2(E,t)\approx \underset{t\to \infty }{\lim }\frac{\Delta X^2(E,t)}{2t}.$$

$$L(E,t)=2\sqrt{\Delta X^2(E,t)},$$

$$\xi (E)=\underset{t\to \infty }{\lim }\frac{L(E,t)}{2\pi }.$$

$$\text{Tr}\left[A\right]\approx ⟨\varphi |A|\varphi ⟩,$$

$$\rho (E)\approx \frac{2}{\Omega }⟨\varphi |\delta (E-H)|\varphi ⟩;$$

$$\rho (E)C_{vv}(E,t)\approx \frac{2}{\Omega }\text{Re}\left[⟨\varphi |U(t)V\delta (E-H)U(t{)}^{†}V|\varphi ⟩\right];$$

$$\rho (E)\Delta X^2(E,t)\approx \frac{2}{\Omega }⟨\varphi |[X,U(t){]}^{†}\delta (E-H)[X,U(t)]|\varphi ⟩.$$

$$U(\pm \Delta t)|\psi ⟩\approx \sum _{m=0}^{N_p-1}(2-{\delta }_{0m})(\mp i{)}^m{J}_m\left(\frac{\widetilde{\Delta t}}{\hslash }\right)T_m(\tilde{H})|\psi ⟩;$$

$$[X,U(\Delta t)]|\psi ⟩\approx \sum _{m=0}^{N_p-1}(2-{\delta }_{0m})(-i{)}^m{J}_m\left(\frac{\widetilde{\Delta t}}{\hslash }\right)[X,T_m(\tilde{H})]|\psi ⟩,$$

$$\tilde{H}=H/\Delta E;$$

$$\widetilde{\Delta t}=\Delta E\Delta t,$$

$$T_m=2\tilde{H}{T}_{m-1}-T_{m-2};$$

$$[X,T_m]=2[X,\tilde{H}]T_{m-1}+2\tilde{H}[X,T_{m-1}]-[X,T_{m-2}];$$

$$T_0=1T_1=\tilde{H};$$

$$[X,T_0]=0,[X,T_1]=[X,\tilde{H}].$$

$$\rho (E)\approx \frac{2}{\pi \Omega \Delta E\sqrt{1-{\tilde{E}}^2}}\sum _{n=0}^{N_m-1}{g}_n(2-{\delta }_{n0})T_n(\tilde{E})C_n^{\mathrm{DOS}};$$

$$\rho (E)C_{vv}(E,t)\approx \frac{2}{\pi \Omega \Delta E\sqrt{1-{\tilde{E}}^2}}\sum _{n=0}^{N_m-1}{g}_n(2-{\delta }_{n0})T_n(\tilde{E})C_n^{\mathrm{VAC}}(t);$$

$$\rho (E)\Delta X^2(E,t)\approx \frac{2}{\pi \Omega \Delta E\sqrt{1-{\tilde{E}}^2}}\sum _{n=0}^{N_m-1}{g}_n(2-{\delta }_{n0})T_n(\tilde{E})C_n^{\mathrm{MSD}}(t).$$

$$C_n^{\mathrm{DOS}}\approx ⟨\varphi |T_n(\tilde{H})|\varphi ⟩;$$

$$C_n^{\mathrm{VAC}}(t)\approx \text{Re}\left[⟨\varphi |U(t)V{T}_n(\tilde{H})U(t{)}^{†}V|\varphi ⟩\right];$$

$$C_n^{\mathrm{MSD}}(t)\approx ⟨\varphi |[X,U(t){]}^{†}{T}_n(\tilde{H})[X,U(t)]|\varphi ⟩.$$

Full text

††journal: Computer Physics Communications

PROGRAM SUMMARY

Manuscript Title: GPUQT: An efficient linear-scaling quantum transport code fully implemented on graphics processing units

Authors: Zheyong Fan, Ville Vierimaa, and Ari Harju

Program Title: GPUQT

Journal Reference:

Catalogue identifier:

Licensing provisions:

*Programming language: * CUDA

Computer: Architectures with CUDA-enabled NVIDIA GPUs with compute capability 2.0 or higher.

Operating system: Linux.

RAM: Needs about 1 to 10 GB device memory and less CPU memory, depending on the size of simulated system.

Number of processors used: One CPU processor and one GPU card.

Keywords: Quantum transport; Linear-scaling; GPU acceleration

Classification: 7.9 Transport properties.

Nature of problem:

Obtain intrinsic electronic transport properties of large systems described by real-space tight-binding Hamiltonians.

Solution method:

The DC conductivity is represented as a time integral of the velocity auto-correlation (VAC) or a time derivative of the mean square displacement (MSD). The calculations achieve linear scaling (with respect to the number of orbitals in the system) computation time and memory usage by using various numerical techniques, including sparse matrix-vector multiplication, random phase approximation of trace, Chebyshev expansion of quantum evolution operator, and kernel polynomial method for quantum resolution operator.

Restrictions:

The number of orbitals is restricted to about 20 million due to the limited amount of device memory in current GPUs.

Running time: About 3 minutes (using a Tesla K40 GPU card) for both examples provided.

1 Introduction

Electrical current can be either viewed as [1] a consequence of an applied electric field as in the classical Boltzmann formalism and the quantum Kubo formalism, or as transmission of charge carriers as in the Landauer-Büttiker formalism. The Landauer-Büttiker formalism, or more generally the non-equilibrium Green’s function formalism, is very versatile and has become the standard method for quantum transport simulations from nanoscale to mesoscale [2, 3, 4, 1, 5]. Combined with the recursive Green’s function technique [6], the Landauer-Büttiker formalism can be used to efficiently simulate relatively narrow systems. An efficient and flexible open source code, Kwant [7], is available for quantum transport simulations based on tight-binding models. As matrix inversion is at the heart of the recursive Green’s function technique, the computational cost generally scales cubically with respect to the width of the system, which severely restricts the application of the method to realistically large 2D and 3D systems. To study large systems, linear-scaling computational cost is desirable. Fortunately, such a linear-scaling method has been developed in the Kubo formalism [8, 9, 10, 11, 12, 13] and has found a lot of applications, especially in quasi-1D and 2D materials [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. For a review, see Ref. [27].

Recently, we have made an efficient GPU implementation [21] of this method. Here, we present our GPU code, which we call GPUQT, and use a few examples to illustrate its usage. GPUQT is fully implemented on the GPU using the CUDA toolkit [28]. Using a single modern graphics card such as Tesla K40, the speedup factor achieved by GPUQT over a serial CPU implementation ranges from one to two orders of magnitude, depending on the problem. Usually, a higher speedup factor can be obtained in a problem with a denser Hamiltonian due to the higher arithmetic intensity, a measure of the amount of floating-point operations relative to the amount of memory accesses required to support those operations. Using GPUQT, one can easily simulate tight-binding systems with millions of sites. To use GPUQT, one has to define a simulation model by specifying the Hamiltonian and current (velocity) operators. The Hamiltonian should be defined in real space and be relatively sparse. As this method is not very suitable for studying ballistic transport properties [21, 14], studied system should also contain one or more types of disorder such that a diffusive regime can be reached within a reasonable computation time.

This paper is organized as follows. We first present the basic theoretical formalisms underlying GPUQT in section 2 and then discuss the numerical techniques crucial for achieving linear scaling in section 3. In section 4, we describe the overall structures of the GPUQT package and specifications in the inputs and outputs. Two examples are presented in section 5 to illustrate the usage of GPUQT. Section 6 summarizes and concludes.

2 Theoretical formalisms

2.1 The Kubo-Greenwood formula for DC conductivity

The Kubo-Greenwood formula [29, 30] for DC diagonal conductivity $\sigma^{\rm KG}(E)$ as a function of the Fermi energy $E$ at zero temperature is

[TABLE]

where $\hbar$ is the reduced Planck constant, $e$ is the electron charge, $\Omega$ is the system volume, $V$ is the velocity operator in the transport direction, $H$ is the Hamiltonian of the system, and Tr denotes the trace. Spin degeneracy is included by the factor 2 in the formula. Linear-scaling evaluation of the Kubo-Greenwood conductivity has been studied early by Thouless and Kirkpatrick [31], later by Mayou[8], and recently by Ferreira and Mucciolo [32].

2.2 Velocity autocorrelation and mean square displacement

In GPUQT, we do not directly calculate the Kubo-Greenwood conductivity $\sigma^{\rm KG}(E)$. Instead, we first calculate one of the two correlation functions, the velocity autocorrelation (VAC) or the mean square displacement (MSD). Both correlation functions are a function of the correlation time $t$.

By Fourier transforming one of the $\delta$ functions in the above formula, one can express the running electrical conductivity $\sigma^{\rm VAC}(E,t)$ as a time integral of the VAC $C_{vv}(E,t)$,

[TABLE]

[TABLE]

[TABLE]

where $V(t)=U^{\dagger}(t)VU(t)=e^{iHt/\hbar}Ve^{-iHt/\hbar}$ is the velocity operator in the Heisenberg representation, and $\rho(E)$ the density of states (DOS). The calculation of the VAC using linear-scaling techniques was first done by Mayou and Khanna [9].

Equivalently, one can expresses the running electrical conductivity as a time-derivative of the mean square displacement (MSD) $\Delta X^{2}(E,t)$,

[TABLE]

[TABLE]

where $X(t)=U^{\dagger}(t)XU(t)$ is the position operator in the Heisenberg representation. This formalism was first proposed by Roche and Mayou [10, 11] and later refined by Triozon et al. [12, 13].

2.3 Transport regimes

The VAC and MSD contain information about the transport regimes. For example, the MSD changes from a quadratic to a linear function of the correlation time during a ballistic-to-diffusive transition and finally saturates, causing the absence of diffusion, if strong (Anderson) localization takes place.

If the transport is diffusive, the VAC usually decays exponentially,

[TABLE]

where $\tau(E)$ is the relaxation time. Then we get the semiclassical conductivity:

[TABLE]

The product of the velocity and the scattering time is the mean free path

[TABLE]

Multiplying this with the velocity gives the diffusivity, $D(E)=v_{x}^{2}(E)\tau(E)$, which, in terms of the MSD, can be understood as an Einstein relation [10, 11]:

[TABLE]

When localization/anti-localization takes place, one usually needs to carefully analyse the behaviour of the running conductivity $\sigma(E,t)$. Quantitative analysis of localization can be facilitated by the use of the MSD. The square root of the MSD,

[TABLE]

serves as a good estimation of the length [17, 18, 21, 22, 23], up to which the electrons propagate. Using this, we can express the running conductivity as a function of the propagating length, $\sigma(E,L)$, rather than the correlation time. In the strongly localized regime, the propogating length will finally saturate and the saturated value is found to be directly related to the localization length [22, 23],

[TABLE]

3 Linear-scaling techniques

Based on the above theoretical formalisms, we see that the quantities that need to be calculated are $\rho(E)$, $\rho(E)C_{vv}(E,t)$, and $\rho(E)\Delta X^{2}(E,t)$. The running conductivity $\sigma(E,t)$ can then be obtained either by a time integral of $\rho(E)C_{vv}(E,t)$, or a time derivative of $\rho(E)\Delta X^{2}(E,t)$. There are a few linear-scaling techniques which work together to achieve linear-scaling computation time and memory usage in the calculations of these quantities.

3.1 Linear-scaling evaluation of the trace

The first approximation is to use a random vector $|\phi\rangle$ to evaluate the trace [33]:

[TABLE]

where $A$ is an arbitrary $N\times N$ matrix operator, and $|\phi\rangle$ is normalized to $N$, $\langle\phi|\phi\rangle=N$. The error introduced by this approximation decreases with increasing $N$, scaling as $\sim 1/\sqrt{N}$ [33]. For a given $N$, the accuracy can also be increased by taking average over independent random vectors. The introduction of the random vector is crucial to achieving linear-scaling, because the major computation will be sparse matrix-vector multiplication, which scales linearly with respect to the system size.

With this approximation, we have

[TABLE]

[TABLE]

[TABLE]

3.2 Linear-scaling evaluation of the time evolution

Both the VAC and MSD formalisms involve a time evolution operator $U(t)$, which is absent from the Kubo-Greenwood formula. After using the random vector approximation, we only need to evaluate the application of the time evolution operator on a vector rather than matrix exponential. The basic idea is to divide the total correlation time into a number of steps. For one time step $\Delta t$ (the time steps need not to be uniform), we have the following Chebyshev polynomial expansions [34, 35, 21]:

[TABLE]

[TABLE]

where $J_{m}$ is the $m$th order Bessel function of the first kind and $T_{m}$ is the $m$th order Chebyshev polynomial of the first kind. Note that $T_{m}$ is defined in the interval $[-1,1]$ and the Hamiltonian and time step have to be scaled in the opposite way:

[TABLE]

[TABLE]

where $\Delta E$ is sufficiently large such that the spectrum of the scaled Hamiltonian $\widetilde{H}$ lies within the interval $[-1,1]$. The Chebyshev polynomial expansions of the time evolution operators can be evaluated up to machine precision and the order of expansion $N_{p}$ needed for achieving this is proportional to the time interval.

The above summations can be efficiently evaluated by using the following recurrence relations ($m\geq 2$) ($T_{m}(\widetilde{H})$ is written as $T_{m}$ for simplicity) [21]:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

3.3 Linear-scaling evaluation of the quantum resolution operator

There are quite a few linear-scaling techniques for approximating the quantum resolution operator $\delta(E-H)$, including the Lanczos recursion method [36, 37], the Fourier transform method [38, 39], and the kernel polynomial method [33]. The Lanczos method is usually not as stable as the kernel polynomial method [33]. In Ref. [21], it has been demonstrated that the Fourier transform method is not as efficient as the kernel polynomial method. We have thus implemented only the kernel polynomial method in GPUQT.

In the kernel polynomial method [33], the quantum resolution operator is approximated by a truncated Chebyshev polynomial expansion, and we can rewrite Eqs. (14-16) as

[TABLE]

[TABLE]

[TABLE]

Here, $C_{n}^{\rm DOS}$, $C_{n}^{\rm VAC}(t)$, and $C_{n}^{\rm MSD}(t)$ are the Chebyshev moments:

[TABLE]

[TABLE]

[TABLE]

A kernel (damping factor) is applied before performing the Chebyshev summation in order to suppress the Gibbs oscillations. For most applications, the Jackson damping [33]

[TABLE]

where $\alpha=1/(N_{m}+1)$, is a good choice. The energy resolution achieved scales as $\delta\sim 1/N_{m}$ [33]. Therefore, to achieve a finer energy resolution, one needs to use a larger $N_{m}$.

4 Using GPUQT

4.1 Compile the code and run the examples

After downloading and unpacking GPUQT, one can see two folders: src and examples. The folder src contains all the source files and a makefile. The folder examples contains two sub-folders with names diffusive and localized, both containing the files make_inputs.cpp and plot_results.m.

To compile the code, simply go to the src folder and type make in the command line. Upon finished, an executable called gpuqt will be created in this folder.

Before running the examples, one has to first go to the diffusive and localized folders and compile (using e.g. g++) and run the make_inputs.cpp code. This will create input files that are needed for running GPUQT.

Then, one needs an extra input file, which we call a “driver input file”, to specify the path(s) of the folder(s) containing the input files. To run the two examples consecutively in a single job, this file should read

examples/diffusive
examples/localized

Suppose that the “driver input file” is named as input.txt and is in the examples folder, one can run the examples using the following command:

src/gpuqt examples/input.txt

4.2 Input files for GPUQT

The input files are used to specify the Hamiltonian of a simulated system and some controlling parameters. All the input files for a simulation should be in a single folder.

In the tight-binding approximation, the Hamiltonian can be written as

[TABLE]

where $H_{mn}$ is the hopping integral between sites (orbitals) $m$ and $n$ and $U_{m}$ is the on-site potential of site $m$. Similarly, the position and velocity operators can be expressed as

[TABLE]

[TABLE]

In the input files, one has to specify four sets of data: 1) the neighbor list structure that determines which hopping integrals are non-vanishing, 2) the non-vanishing hopping integrals, 3) the on-site potentials, and 4) the positions $X_{m}$ of the sites projected onto the transport direction.

4.2.1 The neighbour.in input file

This file specifies the topology of the problem using a neighbour list. This will be used to build the sparse Hamiltonian. The first line should have two integer numbers. The first number is the total number of sites in the simulated system. The second number is the maximum possible number of nonzero hopping integrals originated from a given site. For example, in a square lattice with nearest-neighbour hopping only, this number can be set as 4. Using a larger number than needed will waste memory. Starting from the second line, the $n$th line contains the number of neighbours and the indices of the neighbouring sites of the $(n-1)$th site.

4.2.2 The hopping.in input file

This is an optional input file, which contains the hopping integrals (the off-diagonal terms in the Hamiltonian). If this file is not prepared, GPUQT assumes that all the hopping integrals between pairs of neighbouring sites (specified in the neighbour.in file) are $-1$. The first line should be either the word real or complex. If the word is real, it means that all the hopping integrals are real numbers. Then, starting from the second line, the $n$th line contains the real hopping integrals between the $(n-1)$th site and its neighbouring sites, and the order should be consistent with that in the neighbour.in file. The file will look like this:

real
real_1 real_2 real_3 ...
...

If the word is complex, it means that not all the hopping integrals are real numbers. Then, each real hopping integral as described above should be substituted by two real numbers, the real and imaginary parts of the complex hopping integral. The file will look like this:

complex
real_1 imag_1 real_2 imag_2 real_3 imag_3 ...
...

The unit of energy is determined by the user. One should consistently use the same unit in other input files such as potential.in, energy.in and para.in.

4.2.3 The potential.in input file

This is an optional input file, which contains the on-site potentials (the diagonal terms in the Hamiltonian). If this file is not prepared, GPUQT assumes that all the on-site potentials are zero. The $n$th line is the on-site potential of the $n$th site.

4.2.4 The position.in input file

This file specifies the coordinates of the sites in the simulated system. This will be used to build the velocity (current) operator. The first line should have two numbers, which are the length of the simulated system in the transport direction and the volume of the system. Be careful with periodic boundary conditions. For example, consider a $1~{}000\times 1~{}000$ regular square lattice with a lattice constant of $a=1$, the length in the $x$ direction should be 1 000, even though the distance between a leftmost site and a rightmost site is only 999. Starting from the second line, the $n$th line is the position component of the $(n-1)$th site in the transport direction. The unit of length is determined by the user. One can either set the lattice constant to 1 or some values in unit of nm or Å. What is important is to be consistent when reporting the results.

4.2.5 The energy.in input file

This file contains the energy points to be considered in the simulations. There is a single column in this file. The first line should be an integer, which is the number of energy points to be read in. Starting from the second line, the $n$th line contains the $(n-1)$th energy value. Note that the method is parallel in energy and using one thousand energy points takes roughly as much time as using a single energy point. Here is an example:

601
-3.00
-2.99
...
0.00
...
3.00

This file tells that there would be 601 energy points to be calculated, from $-3$ to $3$, with a spacing of $0.01$.

4.2.6 The time$\textunderscore$step.in input file

This file contains the time steps in the VAC and/or MSD calculations. There is a single column in this fie. The first line should be an integer, which is the number of time steps to be read in. Starting from the second line, the $n$th line is the $(n-1)$th time step. Here is an example:

20
1
2
...
19
20

This file tells that there would be 20 (non-uniform) time steps, from $t_{0}$ to 20 $t_{0}$, with a spacing of $t_{0}$. One should note that the data here are the time steps, not the cumulative times. The cumulative times for this example should be $t_{0}$, $3t_{0}$, $6t_{0}$, $10t_{0}$, $\cdots$. The unit of time, $t_{0}$ is fixed by the energy unit, as we set the reduced Planck constant to 1 in GPUQT. Suppose the unit of energy is $\gamma$, the unit of time is then $t_{0}=\hbar/\gamma$.

4.2.7 The para.in input file

This file contains some additional parameters to define the simulation. In this input file, blank lines are ignored. Each non-empty line starts with a keyword possibly followed by one or more parameters. The valid keywords and their parameters are

calculate_vac

This keyword does not need any parameter. If this keyword appears, the VAC will be calculated. Otherwise, the VAC will not be calculated. 2. 2.

calculate_msd

This keyword does not need any parameter. If this keyword appears, the MSD will be calculated. Otherwise, the MSD will not be calculated. If both this the above keywords are absent, there is no need to prepare the time_step.in input file. 3. 3.

number_of_random_vectors $N_{r}$

This keyword needs one parameter, which is the number of random vectors $N_{r}$ used in the simulation. If this keyword is absent, the default value $N_{r}=1$ will be used. If you want to use 10 random vectors for a given problem, you can either set this number to 10, or set it to 1 and then run the simulation 10 times. Increasing $N_{r}$ can improve the accuracy of the results. 4. 4.

number_of_moments $N_{m}$

This keyword needs one parameter, which is the number of Chebyshev moments $N_{m}$ used in the kernel polynomial method. If this keyword is absent, the default value $N_{m}=1000$ will be used. A larger $N_{m}$ gives a finer energy resolution and one usually needs to test the effects of this parameter. 5. 5.

energy_max $\Delta E$

This keyword needs one parameter, which is a scaling parameter $\Delta E$ used to scale the Hamiltonian. The scaled Hamiltonian $H/\Delta E$ must have all of its eigenvalues lying within the interval $[-1,1]$. If this keyword is absent, the default value $\Delta E=10$ will be used. Using a value larger than needed will only effectively reduce the energy resolution, but using a value smaller than needed will cause big problems as this will lead to calculating the square roots of negative numbers.

4.3 Output files of GPUQT

We now describe the data format of the output files produced by running GPUQT. We note that for all the output files, results from a new simulation will append to, rather than overwrite existing data.

4.3.1 The dos.out file

The $n$th column of this file corresponds to the value of $\rho(E_{n})$ at the $n$th energy point $E_{n}$ specified in the energy.in file. Each row corresponds to the results obtained by using one random vector. The unit of DOS is $1/\gamma/a^{2}$ in 2D and $1/\gamma/a^{3}$ in 3D, where $\gamma$ is the unit of energy and $a$ is the unit of length.

4.3.2 The vac.out file

The $n$th column of this file corresponds to the value of $\rho(E_{n})C_{vv}(E_{n},t)$ at the $n$th energy point $E_{n}$ specified in the energy.in file. If the number of time steps specified in the time_step.in file is $N_{t}$, the first $N_{t}$ rows correspond to the results obtained by using one random vector. Integrating this quantity with respect to time gives the running electrical conductivity. In 2D, the unit of conductivity is $e^{2}/\hbar$. In 3D, the unit is $e^{2}/\hbar/a$, where $a$ is the unit of length. As expected, the unit in 3D can be converted to S/cm.

4.3.3 The msd.out file

The $n$th column of this file corresponds to the value of $\rho(E_{n})\Delta X^{2}(E_{n},t)$ at the $n$th energy point $E_{n}$ specified in the energy.in file. If the number of time steps specified in the time_step.in file is $N_{t}$, the first $N_{t}$ rows correspond to the results obtained by using one random vector. Taking derivative of this quantity with respect to time and then dividing by 2 gives the running electrical conductivity.

5 Examples

In this section, we present two examples to illustrate the usage of GPUQT. Although this method has been mostly used to study graphene-based materials [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26], any system with an appropriate real-space tight-binding Hamiltonian can be treated. Here, for pedagogical purposes, we consider the Anderson model of square lattice. The tight-binding Hamiltonian reads

[TABLE]

where $\langle m,n\rangle$ means a pair of sites which are nearest neighbors of each other. The on-site potentials $U_{m}$ take values uniformly distributed in $[-W/2,W/2]$, where $W$ is the strength of the Anderson disorder.

5.1 The diffusive regime in a 2D lattice

In this example, we consider a $2000\times 2000$ square lattice with $W=4\gamma$, using periodic boundary conditions in both directions. This generally represents an effectively 2D system, although one usually needs to check possible finite-size effects. The input files can be generated by compiling and running the code make_inputs.cpp prepared in the examples/diffusive folder. Using $\gamma$ as the energy unit, all the nonzero hopping integrals have the value of $-1$. Therefore, the hopping.in file is not needed. A uniform time step of $0.2\hbar/\gamma$ is used and the number of time steps is 20, which gives a total correlation time of $4\hbar/\gamma$. The energy points considered range from $-6\gamma$ to $6\gamma$, with an interval of $0.02\gamma$. The para.in file reads:

energy_max 6.1
calculate_vac
calculate_msd

which means that $N_{m}=1000$ (default value), $\Delta E=6.1$, $N_{r}=1$ (default value), and both VAC and MSD will be calculated.

After running this example, three output files, dos.out, vac.out, and msd.out will be generated. Running the MATLAB script plot_results.m prepared in the same folder will produce the graphs shown in Fig. 1. The DOS in Fig. 1 (a) shows that the van Hove singularity at the band center in ordered 2D square lattice disappears in the presence of Anderson disorder. With increasing correlation time $t$, the VAC [Fig. 1 (b)] decays, while the MSD [Fig. 1 (c)] changes from a quadratic to a linear function, both indicating a ballistic-to-diffusive transition. The running conductivities calculated from the VAC and the MSD, $\sigma^{\rm VAC}(E,t)$ and $\sigma^{\rm MSD}(E,t)$, are equivalent to each other, as demonstrated in Fig. 1 (d) for the band center $E=0$. Up to a correlation time of $4\hbar/\gamma$, diffusive transport has been achieved for all the energy points, as can be seen from Fig. 1 (e). Further increasing the correlation time (propagating length) will bring the transport into the localized regime, leading to decreasing running conductivity. Therefore, it is reasonable to take the maximum value (there are exceptions [23], though) at each energy $E$ as the semiclassical conductivity $\sigma_{sc}(E,t)$, which is shown in Fig. 1 (f). From the semiclassical conductivity and the electron group velocity, which is simply square root of the VAC at zero correlation time, one can obtain the relaxation time and the mean free path using Eqs. (8) and (9).

5.2 The localized regime in a 1D chain

In this example, a 1D disordered chain of length $4~{}000~{}000$ and with $W=4\gamma$ is considered, using periodic boundary conditions. The input files can be generated by compiling and running the code make_inputs.cpp prepared in the examples/localized folder. Again, all the nonzero hopping integrals have the value of $-1$ in unit of $\gamma$ and the hopping.in file is not needed. The energy points considered range from $-4\gamma$ to $4\gamma$, with an interval of $0.02\gamma$. As explained in Ref. [21], the VAC formalism is not as practical as the MSD formalism in the localized regime. Also, the time steps do not need to be uniform. The time_step.in for this example reads:

30
0.1
0.2
...
1.0
1
2
...
10
10
20
...
100

That is, we use larger and larger time steps, which can capture both the ballistic-to-diffusive transition and the diffusive-to-localized transition. The para.in file reads:

number_of_random_vectors 2
number_of_moments        500
energy_max               4.1
calculate_msd

which means that $N_{m}=500$, $\Delta E=4.1$, $N_{r}=2$, and only MSD will be calculated.

After obtaining the output files by running GPUQT, one can run the MATLAB script plot_results.m prepared in the same folder to get the graphs shown in Fig. 2. The DOS in Fig. 2 (a) shows smoothed Hove singularities at the band edges in 1D chain. The running conductivity in Fig. 2 (b) clearly shows a maximum at each energy corresponding to $\sigma_{sc}(E,t)$ and the vanishing of conductivity at large correlation time due to Anderson localization. As discussed in Ref. [22], the conductance $g(E,L)=\sigma(E,L)A/L$ ($A$ is the width of the system, which is 1 here) decays exponentially in an appropriate range of the propagating length $L(E)$:

[TABLE]

as shown in Fig. 2(c). The localization length $\xi(E)$ is equivalent to that defined in Eq. (12), as can be seen from Fig. 2(d).

6 Summary and Conclusions

We have presented GPUQT, an efficient CUDA code suitable for studying intrinsic quantum transport properties of large systems described by real-space tight-binding Hamiltonians. Although we have used 2D square lattice and 1D chain with Anderson disorder to illustrate the usage, the inputs to GPUQT are made as general as possible such that many more realistic problems can be studied. The current version is only able to calculate the diagonal conductivity. Hall conductivity [40] and spin relaxation time [41, 42] can also be calculated within the same framework and their implementation will be included in a future version. The code and its updating can be accessed from GitHub [43].