Tutorial 6 Transport with Siesta

Introduction

This tutorial will show you how to setup the Hamiltonian matrix H and the overlap matrix S needed to perform and electron transport calculation.

Important !

If you do not already have the file iosiesta.so in ASE/Transport/Siesta, you will need to create it. The file contains a library which construct the Hamiltonian and overlap matrix from a Siesta calculation and makes them available as arrays in python.

You need to have numpy (www.numpy.org) installed. Go to ASE/Transport/Siesta/f_extension and run:

f2py -c --help-fcompiler

to get a list of available and supported Fortran compilers. Then run the following with "compilervendor" substituted by one of the available compilers:

f2py -c io.f90 hs.f90 -m iosiesta --fcompiler=compilervendor

and finally copy the created library iosiesta.io to ASE/Transport/Siesta/:

cp iosiesta.so ../

CO adsorbed on a 1d Au chain

Lead Calculation

Au wire DFT calculation:

from ASE.Calculators.SIESTA import SIESTA
from ASE import Atom, ListOfAtoms
atoms = ListOfAtoms([Atom('Au',(0.0,0.0,0.0))],cell=(2.9,10.0,10.0),periodic=True)
calc = SIESTA(out='out_Au_wire',xc='PBE',kpts=[65,1,1])
calc.SetMeshCutoff(100)
calc.SetSaveHS(True)
atoms.SetCalculator(calc)
atoms.GetPotentialEnergy()

Setting up the H and S matrices for the Au wire:

from ASE.Transport.Siesta.setup_hs import setup_hs_lead
hs_lead = setup_hs_lead(jobname='out_Au_wire',d=0,gamma=False)
hs_lead.Write(filename='HS_lead.nc',repeat=9,nkpt=51,skpt_c=(0.0,0.0,0.0),spin=1)

Here we setup a H and S matrix for a system corresponding to having repeated the supercell in the DFT calculation 9 times (repeat=9) in the x-direction (d=0), e.g d=2 corresponds to the z-direction. The nkpt argument is the number of k point sampling the transport direction Brillouin zone, and should be an odd integer. skpt_c is a transverse k-point in scaled coordinates, e.g. skpt_c is in the range (-0.5..0.5,-0.5..0.5,-0.5..0.5). If a spin-polarized calculation is performed spin=1 and spin=2 can be used to distingiush between the spin channels.

Central region calculation

Au wire with CO ontop (atomic positions take from PRL 93, 096404 (2004)):

au5_co_loa

The atoms for the scattering region should include two principal layers (of the lead). The first principal layer should be the first atoms in the ListOfAtoms and the second principal layer should be the last atoms. Furthermore it is important (for now), that the order of the position of the atoms within the principal layers are the same as in the lead. Script:

from ASE.Calculators.SIESTA import SIESTA
from ASE import Atom, ListOfAtoms
a = 2.9
b = 5.0
atoms = ListOfAtoms([Atom('Au',(0*a,b,b)),
                     Atom('Au',(1*a,b,b)),
                     Atom('Au',(2*a,b,b)),
                     Atom('Au',(3*a,b,b)),
                     Atom('Au',(4*a,b,b+0.2)),
                     Atom('C' ,(4*a,b,b+1.96+0.2)),
                     Atom('O' ,(4*a,b,b+1.96+1.15+0.2)),
                     Atom('Au',(5*a,b,b)),
                     Atom('Au',(6*a,b,b)),
                     Atom('Au',(7*a,b,b)),
                     Atom('Au',(8*a,b,b))],
                    periodic=True,
                    cell=(9*a,2*b,2*b))
calc = SIESTA(out='out_Au_CO',xc='PBE',kpts=[2,1,1])
calc.SetMeshCutoff(100)
calc.SetNumberPulay(10)
calc.SetSaveHS(True)
atoms.SetCalculator(calc)
atoms.GetPotentialEnergy()

Setting up the H and S matrices for the central region can be done like this:

from ASE.Transport.Siesta.setup_hs import setup_hs_scat

prin = 3 * 15
hs_scat = setup_hs_scat(jobname='out_Au_CO',prin=prin,d=0,basis_size='DZP')
hs_scat.Write(filename='HS_scat.nc',leadfilename='HS_lead.nc',spin=1)

where prin is number of basis functions within a principal layer and should be chosen to give nearest neighbour principal layer interactions only. Here we use a DZP basis set, which results in 15 basis functions (bf) per Au atom. The size of the principial then corresponds to 3 atoms, i.e 3atom * 15 bf/atom = 45 bf.

Transport Calculation

Transport script:

import Numeric as num
from ASE.Transport.ScatteringHamiltonianMatrix \
     import ScatteringHamiltonianMatrix
from ASE.Transport.LeadHamiltonianMatrix import LeftLeadHamiltonianMatrix,\
                                                RightLeadHamiltonianMatrix
from ASE.Transport.MatrixIO import ReadHSMatrix
from ASE.Transport.TransmissionTools import TransmissionTool

prin = 3 * 15
h,s,k = ReadHSMatrix('HS_scat.nc')
h2,s2,k2 = ReadHSMatrix('HS_lead.nc')

hmat = ScatteringHamiltonianMatrix(leftprincipallayer=prin,
                                   rightprincipallayer=prin,
                                   hamiltonianmatrix=h,
                                   overlapmatrix=s)


left = LeftLeadHamiltonianMatrix(principallayer=prin,
                                 hamiltonianmatrix=h2,
                                 overlapmatrix=s2)

right = RightLeadHamiltonianMatrix(principallayer=prin,
                                   hamiltonianmatrix=h2,
                                   overlapmatrix=s2)

energies = num.arange(-4.0,3.0,0.02)
trans = TransmissionTool(hmat,left,right,energies,prin)
trans.AlignFermiLevels()
trans.InitScatteringGreenFunction()
trans.GetLeftLeadGreenFunction().SetInfinitesimal(1.0e-3)
trans.GetRightLeadGreenFunction().SetInfinitesimal(1.0e-3)
trans.GetScatteringGreenFunction().SetInfinitesimal(1.0e-4)
trans.UpdateTransmission()
trans.UpdateDOS()
trans.UpdateEigenChannels()
trans.UpdateSpectralFunctions()
trans.UpdateToNetCDFFile('trans.nc')

Bipyridine 1,4 dithiolate between Au(111) leads

In this section we show how to perform af calculation for a system with more realistic leads. We will consider the bipyridine 1,4 dithiolate sandwiched between Au(111) leads.

Lead calculation

The python script below will peform: i) DFT calculation for a Au(111) lead with 3x3 atoms in the supercell surface plane. ii) Construction of the H and S matrices for a transverse k-points corresponding to 4x4 Monkhorstpack sampling (i.e. 8 irreducible k-points):

from ASE.IO.xyz import ReadXYZ
from ASE.Calculators.SIESTA import SIESTA

atoms = ReadXYZ('Au111_3x3.xyz')
atoms.SetBoundaryConditons((1,1,1))
atoms.SetUnitCell([[21.719917, 0.000000, 0.000000],
                   [ 0.000000, 8.867119, 0.000000],
                   [ 0.000000, 4.433560, 7.679150],fix=True)

calc = SIESTA(out='out_Au111',xc='PBE',kpts=[2,6,6])
calc.SetMixingWeight(0.05)
calc.SetNumberPulay(7)
calc.SetMeshCutoff(200)
calc.SetFDFEntry('PAO.BasisSize','SZP')
calc.SetFDFEntry('PAO.EnergyShift','0.01 Ry')
calc.SetSaveHS(True)
atoms.SetCalculator(calc)
atoms.GetPotentialEnergy()

#Write the lead Hamiltonian and overlap matrix
from ASE.Transport.Siesta.io import io
from ASE.Transport.MatrixIO import WriteHSMatrix

ios = io('out_Au111')
ios.ReadHS()
ios.ReadLOA()
Ef = ios.GetFermiLevel()

skpt_kc = [(0.0,0.375, 0.375),
           (0.0,0.375, 0.125),
           (0.0,0.375,-0.125),
           (0.0,0.375,-0.375),
           (0.0,0.125, 0.375),
           (0.0,0.125, 0.125),
           (0.0,0.125,-0.125),
           (0.0,0.125,-0.375)]

for i,skpt_c in enumerate(skpt_kc):
    h_nn,s_nn = ios.GetHSMatrix(skpt_c=skpt_c)
    h_nn -= Ef * s_nn
    WriteHSMatrix('HS_lead%i.nc' % i,h_nn,s_nn,skpt_c)

Au111_3x3.xyz

Central region calculation

For the central region, we consider a BDT molecule with the sulfur atoms at the fcc hollow sites as shown in the figure below.

Au111_3x3_BDT_geom

Central region DFT calculation and construction of the Hamiltoninan and overlap matrices:

from ASE.Calculators.SIESTA import SIESTA
from ASE.IO.xyz import ReadXYZ

atoms = ReadXYZ('Au111_3x3_BDT.xyz')
atoms.SetBoundaryConditions((1,1,1))
atoms.SetUnitCell([[24.188090, 0.000000, 0.000000],
                   [ 0.000000, 8.867119, 0.000000],
                   [ 0.000000, 4.433560, 7.679150],fix=True)


calc = SIESTA(out='out_Au111_BDT',xc='PBE',kpts=[1,4,4])
calc.SetMixingWeight(0.05)
calc.SetNumberPulay(7)
calc.SetMeshCutoff(200)
calc.SetFDFEntry('PAO.BasisSize','SZP')
calc.SetFDFEntry('PAO.EnergyShift','0.01 Ry')
calc.SetSaveHS(True)
atoms.SetCalculator(calc)
atoms.GetPotentialEnergy()

#Write out the H and S matrices
from ASE.Transport.Siesta.setup_hs import setup_hs_scat
prin = 3 * 4 * 9
hs_scat = setup_hs_scat(jobname='out_Au111_BDT',prin=prin,d=0,
                        basis_size='SZP')

for i in range(8):
    hs_scat.Write(filename='HS_scat%i.nc' % i,
                  leadfilename='HS_lead%i.nc' % i,
                  spin=1)

Au111_3x3_BDT.xyz

Transport calculation

The transmission function should be averaged over transverse k-points. A single transverse k-point transport calculation can done like this:

import Numeric
from ASE.Transport.ScatteringHamiltonianMatrix \
     import ScatteringHamiltonianMatrix
from ASE.Transport.LeadHamiltonianMatrix import LeftLeadHamiltonianMatrix,\
                                                RightLeadHamiltonianMatrix
from ASE.Transport.MatrixIO import ReadHSMatrix
from ASE.Transport.TransmissionTools import TransmissionTool
num = Numeric

prin = 3 * 4 * 9
h,s,k = ReadHSMatrix('HS_scat0.nc')
h2,s2,k2 = ReadHSMatrix('HS_lead0.nc')

hmat = ScatteringHamiltonianMatrix(leftprincipallayer=prin,
                                   rightprincipallayer=prin,
                                   hamiltonianmatrix=h,
                                   overlapmatrix=s)


left = LeftLeadHamiltonianMatrix(principallayer=prin,
                                 hamiltonianmatrix=h2,
                                 overlapmatrix=s2)

right = RightLeadHamiltonianMatrix(principallayer=prin,
                                   hamiltonianmatrix=h2,
                                   overlapmatrix=s2)

energies = num.arange(-10.0,4.0,0.02)
trans = TransmissionTool(hmat,left,right,energies,prin)
trans.AlignFermiLevels()
trans.InitScatteringGreenFunction()
trans.GetLeftLeadGreenFunction().SetInfinitesimal(1.0e-1)
trans.GetRightLeadGreenFunction().SetInfinitesimal(1.0e-1)
trans.GetScatteringGreenFunction().SetInfinitesimal(1.0e-3)
trans.UpdateTransmission()
trans.UpdateDOS()
trans.UpdateEigenChannels()
trans.UpdateSpectralFunctions()
trans.UpdateToNetCDFFile('trans0.nc')

The results of the transmission calculation is shown in the figure below. The dashed dotted (red) line corresponds to a single k-point: (0.0,0.375, 0.375), while the full line is result of averaging over the 8 irreducible k-points of a 4x4 Monkhorstpack grid.

Au111_3x3_BDT_T

xmgrace file:

Au111_3x3_BDT_T.agr

ase2: ASE_Tutorial_6 (last edited 2010-10-20 09:11:16 by localhost)