

In [ ]:

    
from __future__ import print_function
import sisl
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Perform a 4-terminal calculation with 2 crossed Carbon chains.
Running a two-terminal calculation with TranSiesta is a breeze compared to running $N>2$-electrode calculations. When performing $N>2$-electrode calculations an endless combination of different applied bias settings become apparent.
This will be reflected in an even more verbose input for TranSiesta to describe all the 4 electrodes, contours, chemical potentials etc.

This example will primarily create the geometries, and then you should perform data analysis.



In [ ]:

    
chain = sisl.Geometry([[0,0,0]], atom=sisl.Atom[6], sc=[1.4, 1.4, 11])
elec_x = chain.tile(4, axis=0).add_vacuum(11 - 1.4, 1)
elec_x.write('ELEC_X.fdf')
elec_y = chain.tile(4, axis=1).add_vacuum(11 - 1.4, 0)
elec_y.write('ELEC_Y.fdf')
chain_x = elec_x.tile(4, axis=0)
chain_y = elec_y.tile(4, axis=1)
chain_x = chain_x.translate(-chain_x.center(what='xyz'))
chain_y = chain_y.translate(-chain_y.center(what='xyz'))



In [ ]:

    
device = chain_x.append(chain_y.translate([0, 0, -chain.cell[2, 2] + 2.1]), 2)
# Correct the y-direction vacuum
device = device.add_vacuum(chain_y.cell[1, 1] - chain_x.cell[1,1], 1)
device = device.translate(device.center(what='cell'))
device.write('DEVICE.fdf')
device.write('DEVICE.xyz')

Exercises

Run the two electrodes (RUN_ELEC_X/Y.fdf).
Run the TranSiesta analyzation step (see fdf flag: TS.Analyze) and determine the optimal pivoting scheme used.
If you are interested you may try to use the worst pivoting scheme and see if it affects the execution time (however this system is very small so the time difference may be very small).
After analyzation and adding the resulting pivoting scheme to RUN.fdf; run the device (RUN.fdf).
Try and extract similar data as done in TB 6. At least plot one of the DOS quantities.
Extend your DOS plot to be orbitally resolved by extracting only subsets of DOS, in this regard also play with the norm keyword, try and plot the DOS per $s$, sum of $p$, etc. for the orbitals on the Carbon atoms.
- A file named siesta.ORB_INDX has been created by Siesta which contains the orbital information per atom, this should give you access to the indices for extraction.



In [ ]:

    
tbt = sisl.get_sile('siesta.TBT.nc')



In [ ]:

Plot bond (vector) currents, below is a skeleton code to do this, look in the sisl manual for extraction of vector current



In [ ]:

    
xy = tbt.geometry.xyz[:, :]
J1 = # fill in the corresponding code here ()
plt.quiver(xy[:, 0], xy[:, 1], J1[:, 0], J1[:, 1]);



In [ ]:

Performing NEGF calculations for $N$-electrodes with an applied bias is extremely difficult and one should first take on this task once the traditional TranSiesta 2-electrode setup is a breeze.

One of the main difficulties in performing good $N$-electrode calculations is the Poisson solution and the boundary conditions. These are more difficult to empose when having more than 2 electrodes or 1 electrode (only 2 electrode setups are easy).

In a 2 terminal device it is obvious how the applied bias is located (a ramp between the two electrodes), however, when dealing with more than 2 elecrodes all electrodes may have a different chemical potential and thus the variations in how to apply the bias becomes infinite.
Your first task is to read through RUN.fdf and figure out which electrode has which chemical potential and draw a small schematic of it.
Calculate the NEGF with a bias of 0.5 eV (please use this command line, for details of options refer to the Siesta manual):
```
siesta -L no_guess_0.5 -V 0.5:eV RUN.fdf > no_guess_0.5.out
```
which applies a bias of 0.5 eV. Read through the output and find the warning which justifies the name no_guess.
There should now be 2 files that lets you plot the bias potential profile, siesta.VH and no_guess_0.5.VH. Use the below method to plot the plane right between the two carbon chains (HINT: calculating the $z$-value at the center between the two chains may be done using a method for the Geometry object)



In [ ]:

    
def plot_grid(grid, plane_dist=1):
    """ Plot a cut through the grid """
    z_index = grid.index(plane_dist, 2)
    x, y = np.mgrid[:grid.shape[0], :grid.shape[1]]
    dcell = grid.dcell
    x, y = x * dcell[0, 0] + y * dcell[1, 0], x * dcell[0, 1] + y * dcell[1, 1]
    fig, ax = plt.subplots(1, 1)
    im = ax.contourf(x, y, grid.grid[:, :, z_index])
    ax.set_xlabel(r'$x$ [Ang]'); ax.set_ylabel(r'$y$ [Ang]')
    ax.set_title('Potential difference [eV]')
    # Also plot the atomic coordinates
    xyz = grid.geometry.xyz
    ax.scatter(xyz[:, 0], xyz[:, 1], 50, 'k', alpha=.6)
    fig.colorbar(im);



In [ ]:

    
# Read in the two different grids:
grid0 = sisl.get_sile('siesta.VH').read_grid()



In [ ]:

    
no_guess = sisl.get_sile('no_guess_0.5.VH').read_grid()
# Specify the geometry so we can add the atoms to the plot
no_guess.set_geometry(device)
plot_grid(no_guess - grid0, 1.) # replace with the correct z-distance

TranSiesta's method of setting boundary conditions for $N$-electrodes is extremely crude since it only fixes the potential on the electrodes. Instead of letting TranSiesta apply the boundary conditions we can provide an external solution to the Poisson problem with proper boundary conditions of the electrodes.
Below is a method to solve the Poisson problem in Python using pyamg. It takes quite some time, so be patient.
You don't need to understand the below script (but I won't hold you back if you want to carefully read it through ;))



In [ ]:

    
G = sisl.Grid
# Define the boundary conditions in the unit-cell
bc = [[G.PERIODIC, G.PERIODIC],
      [G.PERIODIC, G.PERIODIC], # even though left/right meet, we can still use periodic because of the fixed boundaries of the electrodes
      [G.NEUMANN, G.NEUMANN]]

def solve(A, b, tol=1e-10):
    """ Don't worry too much about this, if you don't have pyamg installed it will fail"""
    import pyamg
    from scipy.sparse.linalg import cg
    print('Initiating solver:')
    ml = pyamg.aggregation.smoothed_aggregation_solver(A, max_levels=1000)
    M = ml.aspreconditioner(cycle='W') # pre-conditioner
    residuals = []
    def callback(x):
        residuals.append(np.linalg.norm(b-A*x))
        print('    itt {} with residual {:.2e} / {:.2e}'.format(len(residuals), residuals[-1], residuals[-1] / residuals[0]))
    x, info = cg(A, b, tol=tol, M=M, callback=callback)
    if info != 0:
        raise ValueError('Could not find solution')
    print('Solved the Poisson problem!!!')
    return x

# Atomic indices for the 4 electrodes
elecs = []
# X-left
elecs.append([1., elec_x, np.arange(elec_x.na)])
# X-right
elecs.append([1., elec_x, np.arange(chain_x.na - elec_x.na, chain_x.na)])
# Y-left
elecs.append([-1., elec_y, np.arange(chain_x.na, chain_x.na + elec_y.na)])
# Y-right
elecs.append([-1., elec_y, np.arange(device.na - elec_y.na, device.na)])

# We use the extra dimension because of idx % shape further down below, i.e. broad-casting
shape = np.array([grid0.shape])
grid = sisl.Grid(shape, bc=bc, geometry=device)

# Now create the pyamg matrix
# Because we already have the "correct" boundary conditions as provided above
# this call will also setup the correct boundary conditions.
A, b = grid.topyamg()

# Find placement
for i, (V, elec, elec_idx) in enumerate(elecs):
    print('Finding indices for electrode {}'.format(i + 1))
    # Note this will also work for *any* skewed electrode cell.
    cube = sisl.Cuboid(elec.sc.cell, center=np.average(device.xyz[elec_idx, :], 0))
    idx = grid.index(cube)
    # fold indices into primary cell
    idx = np.unique(idx % shape, axis=0)
    indices = grid.pyamg_index(idx)
    grid.pyamg_fix(A, b, indices, V)

# Ensure we have memory enough (for really large systems
# the sparse matrix may be *VERY* big).
# small trick to not have the solution array twice.
del grid.grid
grid.grid = solve(A, b).reshape(shape.ravel())
grid.write('V.TSV.nc')

After having solved the Poisson problem an additional file is created V.TSV.nc. It contains the Poisson solution with the appropriate boundary conditions. Now run TranSiesta like this:
```
siesta -L guess_0.5 -V 0.5:eV -fdf TS.Poisson:V.TSV.nc RUN.fdf > guess_0.5.out
```
Once completed there will be two solutions of the same applied bias no_guess_0.5 and guess_0.5. It is instructive to compare the two:
1. Compare the output of the two runs, is there a difference in the SCF?
2. Plot the potential cut for the guess_0.5 file as done above for no_guess_0.5, how do they differ? Try and plot a few other plane-cuts.
3. Calculate the transmission for both runs using these commands:
```
tbtrans -L guess_0.5 -V 0.5:eV -D guess_0.5 RUN.fdf > TBT_guess_0.5.out
tbtrans -L no_guess_0.5 -V 0.5:eV -D no_guess_0.5 RUN.fdf > TBT_no_guess_0.5.out
```
  Plot the two transmissions and density of states, are they the same?
If you have VMD/XCrySDen installed try out the command-line utility sgrid to create a cube file containing the potential profile, this page may be useful.



In [ ]: