Problem in transport calculations with delta Hamiltonian files

[nils-wittemeier]

I am a master student and work on electronic transport through single and double wall carbon nanotubes. I have successfully used sisl to create and save the Hamiltonian for electrodes and devices, and run both bias and zero-bias calculations. Now I would like to add a Gate, as a homogeneous electric field, to a calculation with small bias (1meV), for this purpose I created delta-Hamiltonian files with the code below. Running TBtrans for all field strengths (with the .fdf file below), yields low currents for low field strengths, this is expected, since the simulated tube (16,0) is semiconducting (see figure attached at the bottom). For bigger field strength, however, the current starts strongly fluctuating.

Is the problem related to the way I create the files with sisl?

import numpy
import sisl

# Reading coordinates from file and creating geometry
sile = sisl.io.get_sile('../TBT/DEVICE.xyz')
geom = sile.read_geometry()

# Extracting x-components of atoms
x = sisl.Hamiltonian(geom)
for ia in geom:
    x[ia, ia] = geom.axyz(ia)[0]

# Creating dH-files for different electric field strengths
E = numpy.linspace(-2,2,1001)
for e in E:
    with sisl.get_sile('dH/E_{:.3f}.dH.nc'.format(e), mode='w') as fh:
        dH = x * e
        fh.write_delta(dH)

# Define the device Hamiltonian
TBT.HS ../../TBT/DEVICE.nc

# Save the self-energy function files
TBT.Elecs.Out-of-core  T
TBT.Elecs.Gf.Reuse     T   #default

TBT.k [1 1 1]
ElectronicTemperature 5 meV

# These are the definitions of the electrodes
%block TBT.Elec.Left
  HS ../../TBT/ELEC.nc
  semi-inf-direction -A3
  electrode-position 1
%endblock TBT.Elec.Left

%block TBT.Elec.Right
  HS ../../TBT/ELEC.nc
  semi-inf-direction +A3
  electrode-position end -1
%endblock TBT.Elec.Right
TBT.Voltage 0 eV

# Chemical potentials
%block TBT.ChemPots
  Left
  Right
%endblock TBT.ChemPots
%block TBT.ChemPot.Left
  mu V/2
%endblock TBT.ChemPot.Left
%block TBT.ChemPot.Right
  mu -V/2
%endblock TBT.ChemPot.Right

# If one is only interested in the IV curves
# one may reduce the integration to the
# energy window defined by the chemical potentials.
# This will greatly increase throughput although the
# transmission will be calculated in a very reduced
# energy range dependent on the bias.
TBT.Contours.Eta 0. eV
%block TBT.Contours
  IV
%endblock
%block TBT.Contour.IV
   from -|V|/2 - 5 kT to |V|/2 + 5 kT
     delta 0.2 kT
      method mid-rule
%endblock

[zerothi]

Your code looks somewhat correct.

1) The field strength is given in units of the system in Ang which may not be what you want? I.e. is x equally distributed on -x to x? Basically, if your first x coordinate is 0 Ang and the farthest atom is at 10 Ang, then the field will go from 0 to 10 Ang * e. This would correspond to a "top-gate" since the largest field-strength is at your "highest" atom. It is difficult to assess whether you do it as you suspect without knowing the x-coordinates of your geometry?

2) The electrodes are not gated, and they probably should be. Otherwise you need to confine the gate in a small region far from the electrodes. Currently you define the gate on all atoms (also electrodes), but the electrode file does not seem to contain the same applied gate? In this case you are effectively attaching pristine electrodes to a gated system at an abrupt interface. This may indeed cause such oscillations because the electrode/device states are constantly aligned/mis-aligned.

3) You don't have a TBT.dH input in your fdf? I guess you do it on the command-line?

[nils-wittemeier]

Thank you for your help.

1. The tube used here is centered around the axis x=y=0, and atoms are equally distributed on -x to x.
2. I assume, in order to gate the electrodes, the input file has to be edited accordingly? Or is there a similar way?
3. Yes I did indeed do it on the command line, like this:

   V="0.001"
   f="E_0.800.dH.nc"
   tbtrans -V "$V:eV" -fdf TBT.dH:../dH/$f ../RUN.fdf````

[zerothi]

1) Good. 2) You need to create a new electrode per dH. So that will be a bit cumbersome. In this case it is actually easier to do separate DEVICE.nc and ELEC.nc files per gate, and it shouldn't take that long. 3) Good.

[nils-wittemeier]

I did what you suggested and created DEVICE.nc and ELEC.nc per gate and the results do look better:

I have, however, two more questions: I had already created the files above for the case without a gate from the output of another code. Thus I read the Hamiltonian from these files. The Hamiltonian returned by sisl.physics.Hamiltonian.read is always non-orthogonal. I assume it is implemented that way because the function is made to read .nc files created by (Trans)SIESTA and the stored Hamiltonian would not be therefore usually not be orthogonal.

1. I there a fast way of creating an new orthogonal Hamiltonian object from a non-orthogonal one with the same components?
2. From the source code I understand, that sisl.physics.Hamiltonian.read takes keyword arguments and passes them on to sisl.io.siesta.siesta_nc.read_hamiltonian and eventually to sisl.io.siesta.sisl.io.siesta.siesta_nc._read_class_spin but they don't seem to be used after. It looks like this was intended to pass arguments, to object instantiation. Is this correct, or is there way to read an orthogonal Hamiltonian?

[zerothi]

Indeed that looks better :)

1) Yes.

Hcsr = H.tocsr(0)
H = H.fromsp(H.geometry, Hcsr)

2) The reason it states it is a non-orthogonal basis is simply that the nc file contains the overlap matrix (identity). So in reality it is an orthogonal basis.

I have just committed a fix that should correctly figure out whether the basis is orthogonal or not. So if you update sisl you should be good to go. :) Let me know how it works!

I will consider this as fixed. If you have other questions, please open a new issue :)

[nils-wittemeier]

I would like to thank you very much for your help.

I tested your fix. It did work.