Conductance of annulene molecule between graphene nanoribbon leads

STEP 1: Relax the coordinates

• opt.py
• zgnr-annulene.xyz

from gpaw import *
from ase import *
from gpaw.transport.calculator import Transport
from ase.io import read


system = read('zgnr-annulene.xyz')                   # reading input file; *.traj format is also accepted
#system = read('zgnr-annulene.traj')

system.cell = [28.0, 11.0, 44.9305]  # assigning system cell size 
system.pbc=(1,1,1)                                   # periodice along x, y and z
#system.center()                                     # centering in the box
#view(system)

# setting up the calculation
calc = GPAW(h=0.2,                                   # real-space grid spacing
               xc='PBE',                             # assigning exchange-correlation type
               basis='dzp',                          # assigning functional type
               kpts=(1,2,1),                        # assigning k-points
               occupations=FermiDirac(0.2, maxiter=2000),  # 1st argument kBT, 2nd argument number of iterations 
               mode='lcao',                          # using LCAO mode
               txt='ni_gnr_ni.txt',                  # output log file
               maxiter=300,                          # number of max scf-iterations
               poissonsolver=PoissonSolver(nn=2),    # using Poisson solver
               mixer=Mixer(0.01, 5, weight=100.0))   # using pulley mixing

# Set calculator
system.set_calculator(calc)

constraint=constraints.FixAtoms(mask=[i<18 or i >271 for i in range(290)])  # constrained to keep atoms fixed
system.set_constraint(constraint)
#qn=BFGS(system,trajectory='ni_gnr_ni.traj')
relax=optimize.BFGS(system,trajectory='zgnr-annulene_new.traj')          # assigning optimization solver type
relax.run(fmax=0.03)                                                     # max force tolerance

STEP 2: Compute zero-bias transmission function

• main.py
• zgnr-annulene.traj

from gpaw import *
from ase import *
from gpaw.transport.calculator import Transport
import numpy as np
from ase.visualize import view


name = 'zgnr'  # Name of the output file

# reading the input geometry in *.traj format; *.xyz format is also accepted
system = io.read('zgnr-annulene.traj')  

nn = np.argsort(system.get_positions()[:,2])
system = system[nn]         # storing coordinates of the geometry to system
R = system.get_positions()  # storing indices of position
Luc = R[36,2]-R[0,2]        # computing length of electrode (which contains 36 atoms) along transport direction

# Extend the scattering region
left = system[:36].copy()   # copying L-electrode to left
right = system[-36:].copy() # copying R-electrode to right
left.translate([0,0,-Luc])  # shifting coordinates to the left
right.translate([0,0,Luc])  # shifting coordinates to the right
C = system.cell             # storing original cell size
C[2,2] += 2*Luc             # extend the cell size
system = left + system + right  # coordinates of the extended cell
system.set_cell(C)

system.set_pbc([1,1,1])  # repeating cell along 3-orthogonal coordinates; [1, 1, 1] means NO repetition
system.center()          # centering the cell into the box

view(system)
write('system.traj',system)  # writing the extended cell 

lead_cell = np.diag(system.cell)  # original length of the electrode
lead_cell[2]=2*Luc                # defining extended electrodes length
pl_atoms1=range(72)               # assigning atoms of the L-electrode
pl_atoms2=range(290,362)          # assigning atoms of the R-electrode 
pl_cell1= lead_cell               # assigning L-elec size
pl_cell2= lead_cell               # assigning R-elec size

basis= 'dz(dzp)'                       # assigning functional to be used


# Defining the parameters/arguments to be used in transport calculation
t = Transport( h=0.2,                               # grid space along x,y and z
               xc='PBE',                            # exchange functional 
               basis=basis,                         # basis set
               kpts=(1, 1, 1),                      # No of k-points for the electrodes
               occupations=FermiDirac(0.1),         # room temperature kBT = 0.025 eV
               mode='lcao',                         # using LCAO mode
               txt=name+'.txt',                     # output file
               poissonsolver=PoissonSolver(nn=2),
               mixer=Mixer(0.1, 5, weight=100.0),   # using pulley mixing
               pl_atoms=[pl_atoms1, pl_atoms2],     # setting up two electrodes
               pl_cells=[pl_cell1, pl_cell2],       # setting up electrodes cells
               pl_kpts=[1, 1, 15],                  # k-points for the scatting region
               edge_atoms=[[0 , 71],[0 , 361]],     # Defining edge atoms in the form [[L_i, R_f][C_i, C_f]]        
               non_sc = True,                       # True = Normal DFT (default), False = NEGF-DFT
               mol_atoms=range(72, 290),            # molecule atoms in the scattering part
               plot_energy_range = [-0.1,0.1],      # min and max energy
               plot_energy_point_num = 201)         # No of energy points
system.set_calculator(t)
t.negf_prepare()
t.non_sc_analysis()                                 # Calculating T(E=0)

STEP 3: Extract transmission function data from the output file obtained in STEP 2

• transmission.py

• On Chimera execute: mpirun gpaw-python transmission.py '0*' > transmission.dat

from gpaw.transport.analysor import Transport_Plotter
import numpy as np
import sys
from pylab import *
if '*' in sys.argv[1]:
    fd=0
    bias_step = int(sys.argv[1].split('*')[0])
else:
    fd=1
    bias_step = int(sys.argv[1])

plotter=Transport_Plotter(fd)
plotter.plot_setup()
tc = plotter.tc(bias_step, s=0)
#tc2 = plotter.tc(bias_step, s=1)
ee=np.linspace(-1,1,201)
#plot(ee, tc, 'b-o')
#plot(ee, tc, 'r-')
#plot(ee, tc2, 'b-')

for n in range(201):
    print ee[n], tc[n] #, tc2[n]