Tutorial 5 Transport

Introduction

This tutorial describes how to use the ASE Transport module for performing realistic calculations of electron transport in nanoscale contacts. The class TransmissionTool allows you to calculate the Landauer Buttiker conductance (more generally: the elastic transmission function) of any system using a Green's function method. The system is described by a Hamiltonian matrix which must be represented in terms of a localized basis set. The classes HoppingParameters and CompositeWannierCells in the ASE.Utilities.Wannier module allows you to construct such Hamiltonians within DFT in terms of a Wannier function basis set.

In the first part of the tutorial, the use of the transport module will be explained. The second part deals with the construction of realistic DFT Hamiltonians in terms of a Wannier function basis.

For more information on the general transport theory see chapter 3 in the textbook: Electronic Transport in Mesoscopic Systems, S. Datta, Cambridge (1995).

For information on the implementation of the Green's function transport scheme as well as the construction of the Hamiltonian in terms of a Wannier function basis set see the paper: Molecular transport calculations with Wannier functions, Thygesen and Jacobsen, Chem. Phys. in press.

Structure of the Hamiltonian matrix

Before one can start a transport calculation a Hamiltonian matrix describing the system is needed. The Hamiltonian must be constructed in terms of a localized basis set such that a division of the system into left lead (L), right lead (R) and a scattering region (S) is possible. It is assumed that there is no coupling between L and R. It is also assumed that the leads are periodic. A principal layer is defined as the smallest number of lead periods such that only principal layers immediately next to each other couple.

The Hamiltonian of the system can be represented by two Hamiltonian matrices: (i) H_scat which describes the scattering region including 1 principal layer on either side (ii) H_L and H_R which describe the left and right leads. These matrices have the generic shape:

         {H_L  V_L^d    0}
H_scat=  {V_L   H_S   V_R}
         {0    V_R^d  H_R} , where d=dagger

H_i=     {H_i  t_i^d}
         {t_i    H_i}, i=L,R

The dimension of H_scat is (N_scat+N_l+N_R)x(N_scat+N_l+N_R).

The dimension of H_L is (2*N_L)x(2*N_L).

The dimension of H_R is (2*N_R)x(2*N_R).

Here N_scat, N_L, N_R are the number of basis functions in S, L, R, respectively.

Often the left and right leads are identical in which case H_L=H_R and N_L=N_R.

One-dimensional TB chain with two levels (Or: simple model of transport through a hydrogen molecule)

We consider the transport of electrons through a model system where the leads are one-dimensional TB chains with next-nearest neighbor coupling and the scattering region consists of two TB sites coupled weakly to the leads. The system can be viewed as a simple model of a hydrogen molecule sandwiched between two electrodes.

First the lead Hamiltonian is constructed (left and right leads are identical). To include second-nearest neighbor coupling the principal layer must be 2 sites (otherwise there will be coupling beyond nearest neighbor principal layers). Thus H_lead must have dimension 4x4:

import Numeric

H_lead=Numeric.zeros([4,4],Numeric.Complex)

# On-site energies are zero
for i in range(4):
    H_lead[i,i]=0.0

# Nearest neighbor hopping is -1.0
for i in range(3):
    H_lead[i,i+1]=-1.0
    H_lead[i+1,i]=-1.0

# Next-nearest neighbor hopping is 0.2
for i in range(2):
    H_lead[i,i+2]=0.2
    H_lead[i+2,i]=0.2

Next, we construct the Hamiltonian describing the scattering region plus 1 principal layer on each side. Since S consists of 2 TB sites (the 2 hydrogen orbitals), and since a principal layer contains 2 TB sites, the dimension of H_scat is 6x6:

H_scat=Numeric.zeros([6,6],Numeric.Complex)

# Principal layers on either side of S
H_scat[[2,|2]]=H_lead[[2,|2]]
H_scat[-2:,-2:]=H_lead[[2,|2]]

# Scattering region
H_scat[2,2]=0.0
H_scat[3,3]=0.0
H_scat[2,3]=0.8
H_scat[3,2]=0.8

# External coupling
H_scat[1,2]=0.2
H_scat[2,1]=0.2
H_scat[3,4]=0.2
H_scat[4,3]=0.2

The two hydrogen orbitals constituting S have onsite energies 0 (like the sites in the leads) and the H-H coupling is 0.8. The coupling between the hydrogen orbitals and the outhermost sites in the leads is 0.2 .

Initialize classes to contain scattering and lead Hamiltonians:

from ASE.Transport.ScatteringHamiltonianMatrix import ScatteringHamiltonianMatrix
from ASE.Transport.LeadHamiltonianMatrix import LeftLeadHamiltonianMatrix,RightLeadHamiltonianMatrix

# Set the principal layer length (no. of sites in a principal layer)
prin=2
hscat=ScatteringHamiltonianMatrix(leftprincipallayer=prin,rightprincipallayer=prin,hamiltonianmatrix=H_scat)
hleft=LeftLeadHamiltonianMatrix(principallayer=prin,hamiltonianmatrix=H_lead)
hright=RightLeadHamiltonianMatrix(principallayer=prin,hamiltonianmatrix=H_lead)

Select the energy grid on which we want the transmission function, and initialize the TransmissionTool class:

energies=Numeric.arange(-3,3,0.02)
trans=TransmissionTools.TransmissionTool(hscat,hleft,hright,energies,prin)
trans.InitGreenFunctions()

Specify which quantities we want to calculate. These include: (1) The transmission function (2) The spectral functions (projected density of states) for orbitals in S (3) The total density of states of S (4) The transmission eigen-channels Here we want (1) and (2):

trans.UpdateTransmission()
trans.UpdateSpectralFunctions()

Finally, perform calculation (the output is stored in the specified file in the NetCDF file format):

t=trans.UpdateToNetCDFFile('TBchain_trans.nc')

Below is shown the resulting transmission function:

1dtrans

There are two peaks at -0.8 and 0.8, respectively, corresponding to the bonding and anti-bonding combinations of the two hydrogen orbitals. Below is plotted the projected DOS (PDOS) for each of the two sites in the scattering region:

1dpdos

The PDOS is, of course, the same for the two hydrogen orbitals. The weight of the PDOS is located on both the bonding and anti-bonding states.

Further analysis of the transmission function

Below it is shown how to apply the methods in the module HamiltonianTools located in ASE.Utilities.Wannier.

First, it is more useful to describe the system in terms of the bonding and anti-bonding states instead of the hydrogen orbitals. To accomplish this we diagonalize the subspace spanned by the two hydrogen orbitals. The two hydrogen orbitals are states number 2 and 3 in the scattering region Hamiltonian (remember that this matrix include two lead sites on either side of the hydrogen molecule):

from ASE.Utilities.Wannier import HamiltonianTools

# Get the scattering Hamiltonian
H_scat=hmat.GetHamiltonianMatrix()
states=[2,3]
# Diagonalize the subspace spanned by states 2 and 3
H2,U,eig=HamiltonianTools.SubDiagonalize(H,states)
print "The renormalized molecular energy levels",eig

H2 will contain the rotated Hamiltonian matrix. U is the unitary matrix implementing the basis change, and eig are the eigenvalues of the Hamiltonian within the subspace. The above lines produces the output:

Indices of new basis functions: [2, 3]
The renormalized molecular energy levels [ 0.8 -0.8]

The first line shows that after the basis change (from hydrogen orbitals to bonding/anti-bonding states) the new basis functions have the same indices as the original ones (this will always be the case). The second line simply shows the energies of the bonding amd anti-bonding states. Not surprisingly they are exactly at the position of the transmission peaks.

Let use see what happens if we perform a transmission calculation with the bonding state removed from the basis set. This can effectively be done by cutting all coupling matrix elements involving this orbital:

H3=HamiltonianTools.CutCoupling(H2,[3])
hmat.SetHamiltonianMatrix(H3)

The last line updates the Hamiltonian matrix.

Now run the calculation again:

trans=TransmissionTools.TransmissionTool(hmat,left,right,en,prin)
trans.InitGreenFunctions()

trans.UpdateTransmission()
trans.UpdateSpectralFunctions()

t=trans.UpdateToNetCDFFile('TBchain_trans.nc')

The resulting transmission function is shown below:

1dtrans2

Note, that the lower peak has disappeared due to the removal of the bonding state. The projected DOS for the bonding and anti-bonding states are shown below:

1dpdos2

The PDOS of the bonding orbital is not seen because it has been decoupled from all other orbitals. The PDOS of the anti-bonding orbital has weight exactly at the position of the transmission peak.

Constructing a Wannier function Hamiltonian

Before reading this section, the reader should consult the Wannier function tutorial.

The example below illustrates how to construct a Hamiltonian matrix in terms of a Wannier function basis. The first part deals with the Hamiltonian of the leads, H_L and H_R, see the section Structure of the Hamiltonian matrix above, while the second part deals with the Hamiltonian of the scattering region, H_scat.

NOTE: The examples given below assume that the output of certain DFT calculations exists and they can therefore not be run without these output files.

NOTE: It is very important that the grid spacing in the DFT calculation for the scattering region is the same as the grid spacing used in the separate DFT calculation for the lead.

Lead Hamiltonian

Since most transport calculations are performed with identical leads we shall only consider this case here.

First step is to perform a DFT calculation for a unit cell describing the lead material. Typically a large k-point grid will be used for this calculation as the lead is a periodic system. To make sure that unoccupied states (to be used for constructing partly occupied Wannier functions) are properly converged, the calculation should be a Harris calculation, i.e. the electron density should be taken as input, and the DFT code should converge all Kohn-Sham eigenvalues.

In the following we refer to the output of the above DFT calculation as out_lead.nc.

Next, a set of Wannier functions for the lead should be constructed. For information on this step see the Wannier function tutorial.

We refer to the output of the Wannier function construction as zibloch_lead.pickle and Wannier_lead. Also we refer to the number of Wannier functions as N_w_lead and the used occupation energy as E_w_lead.

We are now ready to construct the lead Hamiltonian. First the Wannier object is initialized:

from ASE.Utilities.Wannier.Wannier import Wannier
atoms = Dacapo.ReadAtoms('out_lead.nc')
calc = atoms.GetCalculator()
wannier = Wannier(numberofwannier=N_w_lead,calculator=calc,occupationenergy=E_w_lead)

wannier.ReadZIBlochMatrix('zibloch_lead.pickle')
wannier.ReadRotation('Wannier_lead')

Next, an instance of the class HoppingParameters is initialized, and the Hamiltonian is written to file:

from ASE.Utilities.Wannier.HoppingParameters import HoppingParameters
hop=HoppingParameters(wannier,couplingradius=7.0)
hop.WriteHSMatrixToNetCDFFile(filename='HS_lead_-0.375x-0.375',supercell=[2,3,3],kpoint=[-0.375,-0.375])

The resulting Hamiltonian will describe a system where the lead unit cell has been repeated 2 times in the transport direction (assumed to be the x-axis) and 3x3 times in the transverse directions. The dimension of this Hamiltonian will be (2*3*3*N_w_lead)x(2*3*3*N_w_lead). The transverse k-point has been chosen to (-0.375,-0.375), in scaled coordinates.

Scattering region Hamiltonian

The first step in setting up a Wannier function Hamiltonian for transport calculation is to perform a DFT calculation for the scattering region. It is important that enough lead material has been included in this calculation for the effective potential to have converged to the bulk (lead) value at the end planes of the supercell. To make sure that unoccupied states (to be used for constructing partly occupied Wannier functions) are properly converged, the calculation should be a Harris calculation, i.e. the electron density should be takes as input, an the DFT code should converge all Kohn-Sham eigenvalues. If transverse k-points are needed, a Harris calculation for each tranverse k-point should be performed (along the transport direction the gamma-point should be used).

In the following we refer to the output of the above DFT calculation as out_scatter_kpt-0.375x-0.375.nc assuming that the Harris calculation has been performed at the k-point (0,-0.375,-0.375), where the first coordinate refers to the transport direction.

The next step is then to construct Wannier functions for the scattering region. For information on this step see the Wannier function tutorial.

In the following we refer to the output of the Wannier function construction as zibloch_scatter.pickle and Wannier_scatter. Also we refer to the number of Wannier functions as N_w_scat and the used occupation energy as E_w_scat.

We are now ready to construct the Hamiltonian H_scat. First the relevant modules are imported, and next the DFT files for the scattering region and the lead are read:

import Numeric
from ASE import Atom, ListOfAtoms
from Dacapo import Dacapo
from ASE.Utilities.Wannier import Wannier,CompositeWannierCells

atoms1 = Dacapo.ReadAtoms("out_scatter_kpt-0.375x-0.375.nc",save_memory=True)
atoms2 = Dacapo.ReadAtoms("out_lead.nc",save_memory=True)

calc1 = atoms1.GetCalculator()
calc2 = atoms2.GetCalculator()

The a Wannier object for each of the two regions is initialized:

wannier1 = Wannier(numberofwannier=N_w_scat,calculator=calc1,occupationenergy=E_w_scat)
wannier1.ReadZIBlochMatrix('zibloch_scatter.pickle')
wannier1.ReadRotation('Wannier_scatter')

wannier2 = Wannier(numberofwannier=N_w_lead,calculator=calc2,occupationenergy=E_w_lead)
wannier2.ReadZIBlochMatrix('zibloch_lead.pickle')
wannier2.ReadRotation('Wannier_lead')

To obtain the required Hamiltonian, H_scat, we must glue one lead principal layer on each side of the scattering region. We now construct three cells: one for each of the leads (these are identical in this case) and one for the scattering region:

cell1=CompositeWannierCells.WannierCell(wannier2,atoms2,repeat=[1,3,3])
cell2=CompositeWannierCells.WannierCell(wannier1,atoms1,repeat=[1,1,1])
cell3=CompositeWannierCells.WannierCell(wannier2,atoms2,repeat=[1,3,3])

The keyword repeat specifies how the unit cell should be repeated. In this case the lead unit cell is repeated once in the transport direction and 3x3 in the transverse plane. Of course, the cells should fit nicely at the interfaces. Finally, we construct an object that holds all three cells, specify the transverse k-point, and run the calculation::

container=CompositeWannierCells.CompositeWannierCell([cell1,cell2,cell3],kpoint=[-0.375,-0.375])
container.WriteHSMatrixToNetCDFFile('HS_-0.375x-0.375.nc')

NOTE: Depending on the size of the system, this calculation should be run in parallel. If strange error messenges are received try to increase the number of processors ;-)

Given HS_-0.375x-0.375.nc and HS_lead_-0.375x-0.375.nc a transport calculation can be performed as described in the first example of this tutorial. The only difference is that the Hamiltonian and overlap matrices must be read from file. This is done as follows (here prin is the number of basis functions in one principal layer):

hscat=ScatteringHamiltonianMatrix(leftprincipallayer=prin,rightprincipallayer=prin,filename='HS_-0.375x-0.375.nc')
hleft=LeftLeadHamiltonianMatrix(prin,filename='HS_lead_-0.375x-0.375')
hright=RightLeadHamiltonianMatrix(prin,filename='HS_lead_-0.375x-0.375.nc')

The transport calculation as well as the analysis part is now completely equivalent to what is described for the simple TB model in the example above.