First TranSiesta bias example.

In this example we will take the system from TS 1 and perform bias calculations on it. Note, however, that applying a bias to a pristine bulk system is non-physical and should thus NEVER be done. TranSiesta will not warn you about this and will happily calculate the non-equilibrium density for any bulk system with an applied bias. This is an extremely important point and once complete with this example you should carefully think through why this is the case.

Bias calculations are very heavy because of a full DFT+NEGF calculation per bias-point.

We will begin with creating the structures.

Exercises

In this exercise you will be familiarized with the input options that define a bias calculation. The input options are extremely elaborate, yet they require little intervention when using default parameters.
As this is your first example of running TranSiesta with applied bias there are a few things you should know:

1. Do not start by performing any $V\neq0$ calculations until you are perfectly sure that your $V=0$ calculation is well converged and well behaved, i.e. a small dQ (see TS 1).
2. When performing bias calculations your are recommended to create a new directory for each bias: TS_<V>.

3. Any bias calculation should be a restart by using the closests bias calculations TranSiesta density matrix. This can be ensured by copying the siesta.TSDE file from the closests bias calculation to the current simulation directory. I.e.

   • First run $V=0$ in TS_0, ensure convergence etc.
   • Second run $V=0.25\,\mathrm{eV}$ in TS_0.25, copy TS_0/siesta.TSDE to TS_0.25/ and start run.
   • Third run $V=0.5\,\mathrm{eV}$ in TS_0.5, copy TS_0.25/siesta.TSDE to TS_0.5/ and start run.
   • etc.
   • $N$th run $V=-0.25\,\mathrm{eV}$ in TS_-0.25, copy TS_0/siesta.TSDE to TS_-0.25/ and start run (note negative bias' can be performed in parallel to positive bias)

4. All the commands required for this example can be executed like this:

    siesta --electrode RUN_ELEC.fdf > ELEC.out
    cd TS_0
    cp ../C.psf .
    siesta ../RUN.fdf > TS.out
    # Check that the charge is converged etc.
    cp siesta.TSDE ../TS_0.5/
    cd ../TS_0.5
    cp ../C.psf .
    siesta -V 0.5:eV ../RUN.fdf > TS.out
    # Check that it has converged...
    cp siesta.TSDE ../TS_1.0/
    cd ../TS_1.0
    cp ../C.psf .
    siesta -V 1:eV ../RUN.fdf > TS.out
    # Check that it has converged...

   After every calculation go through the output to ensure everything is well behaved. Note that the output of a bias calculation is different from a non-bias calculation, they are more detailed.

5. An additional analysis (before going to the transport calculations) is to calculate the potential drop in the junction. In sisl this is easy:

   v0 = sisl.Grid.read('TS_0/ElectrostaticPotential.grid.nc')
   vd = (sisl.Grid.read('TS_0.5/ElectrostaticPotential.grid.nc') - v0)

   vd then contains the potential profile (in eV). To save it as a linear average bias file (remember transport is along first lattice vector) you can execute the following:

   vd = vd.average(1).average(2)
   dv = (vd.dcell[0, :] ** 2).sum() ** .5
   sisl.io.tableSile('potential_0.5.dat', 'w').write_data(dv * np.arange(vd.shape[0]), vd.grid[:, 0, 0])

This completes all non-equilibrium calculations for this example. However, we have only calculated the non-equilibrium density and thus, the non-equilibrium Hamiltonian. We still need to calculate the transport properties for all bias'. Basically we can only calculate the transport properties at the calculated bias values, but generally we are interested in a full $I(V)$ curve.
As a user, one has three options:

1. Calculate $I(V)$ for the calculated biases $V$ and perform an interpolation of $I(V)$, or
2. Interpolate the Hamiltonian to calculate $I(V)$ for all the required biases, or
3. Calculate all non-equilibrium Hamiltonians.

The first option is by far the fastests and easiest with a sometimes poor accuracy; the second option is relatively fast, and drastically improves the accuracy; while the last option is the most accurate but may sometimes be non-feasible due to insufficient computational resources.

In the following we will calculate all transmissions using option 2. Look in the manual for the options regarding the interpolation (there are two interpolation methods).
Go through RUN.fdf and find the respective block that tells TBtrans to interpolate the Hamiltonian, also notice how the energy-grid is defined in TBtrans. You will notice that this is the fastest way to calculate the $I(V)$ curve for any bias, it however, will not calculate any physical quantities outside the bias window.

Now complete the exercise by running TBtrans for $V\in\{0, 0.1, \dots, 1\}$ eV. Note that instead of changing the applied bias in the fdf-file, one can do:

tbtrans -V 0.4:eV RUN.fdf

to apply $V=0.4\,\mathrm{eV}$, any fdf-flag specified on the command line has precedence! The : is to denote an effective space, otherwise you will have to encapsulate in quotation marks tbtrans -V "0.4 eV" RUN.fdf.

If you do not want to run the commands manually, you may use this loop command:

for V in $(seq 0 0.1 1) ; do
   d=TBT_${V//,/.}
   mkdir $d
   cd $d
   tbtrans -V "${V//,/.}:eV" ../RUN.fdf > TBT.out
   cd ../
done

TIME: A remark on this exercise. Think why applying a bias to a bulk system is wrong. If you can't immediately figure this out, try and create a longer system by replacing device = elec.tile(3, axis=0) with, say: device = elec.tile(6, axis=0) and redo the calculation for a given bias. Then compare the potential profiles.

Plot the transmissions

Calculate the current for all $V$, then plot it.

Why is the current $0$ for $V<0.8$?
Hint: see TB 1.