 |
 |
|
| Pi model for carbon | Bibliography |
Table of Contents
Next, it is instructive to consider a model with more than one element, to see how that changes things. A commonly used Slater-Koster parametrization for studying Si nanowires is the sp3d5s* model by Zheng et al. [2] which extends an earlier parametrization of Si and Ge [3], which in turn builds on the original sp3d5s* models by Jancu et al. [4]. (The latter is available as a built-in basis set in ATK.)
It is useful to start by defining the model for Si, since it is a little bit more involved than the simple pi-model of the previous chapter, and also quite useful on its own. Then the hydrogen matrix elements will be added to the model; these have a very specific purpose in this model, namely to remove spurious states due to dangling Si bonds on the surface of the nanowire. So, they actually just constitute a numerical trick, as discussed in detail in Ref. [2].
As before, the onsite and offsite matrix elements need to be defined, and the main question is how the latter scale as a function of the distance.
Starting with the onsite parameters, the ionization potentials can be
taken directly from Ref. [3]. There
are two occupied orbitals (3s and 3p), and two unoccupied ones
(3d and an excited s-type orbital, commonly denoted s* in the literature).
Since spin will not be considered (to be precise, the model does contain
spin-orbit interaction, but this is not yet implemented in ATK),
the spin split parameter can be left out. As before, the onsite Hartree
shifts can be extracted from the ATK 
database; for the unoccupied
orbitals 3d and s* will be use values from occupied the valence orbitals:
Es = -2.15168*eV
Ep = 4.22925*eV
Ed = 13.78950*eV
Es1 = 19.11650*eV
onsite = SlaterKosterOnsiteParameters(
element = Silicon,
angular_momenta = [0,1,2,0],
occupations = [2,2,0,0],
ionization_potential = [Es,Ep,Ed,Es1],
onsite_hartree_shift = ATK_U(Silicon,['3s','3p','3p','3s'])
)
This is a nearest-neighbor model, but the scaling function for the matrix elements is fitted to allow for at least somewhat irregular structures. It is commonly assumed that the matrix elements ijk(d) follow a power law
where 
is the ideal neighbor distance.
Harrison [5] originally used a general
exponent of 
for all matrix elements 
;
later work, like Ref. [4],
extends this to fit an individual exponent for each matrix element.
The scaling function can be tabulated in an interval around the nearest-neighbor distance using the following recipe:
# Nearest and second nearest neighbor distances # where a is the cubic lattice constant d_n1 = sqrt(3.0)/4.0*a d_n2 = sqrt(2.0)/2.0*a # Interpolate distances and energies epsilon = linspace(-0.10, 0.10, 21) distances = [ d_n1*(1.0+x) for x in epsilon ] + [ 0.5*(d_n1+d_n2) ] ijk = [ ijk0/(1.0+x)**n_ijk for x in epsilon ] + [ 0.0*eV ]
Here ijk0 is the matrix element from the article,
and the cut-off distance is taken at the midpoint between the first
and second nearest neighbors (the separate element added to the list)
where the corresponding matrix element is set to zero.
This way a power law scaling is obtained around the nearest-neighbor
distance d_n1.
Now, you just need to collect the parameters from the article:
si_si_sss = [ -1.95933*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV] si_si_s1s1s = [ -4.24135*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV] si_si_ss1s = [ -1.52230*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV] si_si_sps = [ 3.02562*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV] si_si_s1ps = [ 3.15565*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV] si_si_sds = [ -2.28485*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV] si_si_s1ds = [ -0.80993*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV] si_si_pps = [ 4.10364*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV] si_si_ppp = [ -1.51801*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV] si_si_pds = [ -1.35554*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV] si_si_pdp = [ 2.38479*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV] si_si_dds = [ -1.68136*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV] si_si_ddp = [ 2.58880*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV] si_si_ddd = [ -1.81400*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
In this code is used for all parameters, but it
is clearly very easy to modify this to have a specific value
for the exponent for each parameter [4].
The final, complete script, with also the Slater-Koster table definition (cf. the rules regarding keyword name in the previous chapter) becomes as follows:
from NanoLanguage import *
'''
Generate a Slater-Koster basis set for silicon based on
T. B. Boykin et al., PRB 69, 115201 (2004)
http://dx.doi.org/10.1103/PhysRevB.69.115201
'''
from math import sqrt
# a is the cubic lattice constant
a = 5.430*Ang
# Nearest and second nearest neighbor distances
d_n1 = sqrt(3.0)/4.0*a
d_n2 = sqrt(2.0)/2.0*a
epsilon = numpy.linspace(-0.10, 0.10, 21)
distances = [ d_n1*(1.0+x) for x in epsilon ] + [ 0.5*(d_n1+d_n2) ]
si_si_sss = [ -1.95933*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_s1s1s = [ -4.24135*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_ss1s = [ -1.52230*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_sps = [ 3.02562*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_s1ps = [ 3.15565*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_sds = [ -2.28485*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_s1ds = [ -0.80993*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_pps = [ 4.10364*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_ppp = [ -1.51801*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_pds = [ -1.35554*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_pdp = [ 2.38479*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_dds = [ -1.68136*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_ddp = [ 2.58880*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_ddd = [ -1.81400*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
Es = -2.15168*eV
Ep = 4.22925*eV
Ed = 13.78950*eV
Es1 = 19.11650*eV
onsite = SlaterKosterOnsiteParameters(
element = Silicon,
angular_momenta = [0,1,2,0],
occupations = [2,2,0,0],
ionization_potential = [Es,Ep,Ed,Es1],
onsite_hartree_shift = ATK_U(Silicon,['3s','3p','3p','3s'])
)
basis_set = SlaterKosterTable(silicon = onsite,
si_si_sss = zip(distances,si_si_sss),
si_si_s1s1s = zip(distances,si_si_s1s1s),
si_si_ss1s = zip(distances,si_si_ss1s),
si_si_sps = zip(distances,si_si_sps),
si_si_s1ps = zip(distances,si_si_s1ps),
si_si_sds = zip(distances,si_si_sds),
si_si_s1ds = zip(distances,si_si_s1ds),
si_si_pps = zip(distances,si_si_pps),
si_si_ppp = zip(distances,si_si_ppp),
si_si_pds = zip(distances,si_si_pds),
si_si_pdp = zip(distances,si_si_pdp),
si_si_dds = zip(distances,si_si_dds),
si_si_ddp = zip(distances,si_si_ddp),
si_si_ddd = zip(distances,si_si_ddd),
)
![[Note]](images/icons/note.png){kind=icon}
|
Note |
|---|---|
|
This basis set contains two s-type orbitals (3s and s*). To distinguish orbitals of the same angular momentum, there are two things to keep in mind:
|
The goal in this chapter is to enable modeling of H-terminated Si nanowires, and now it is time to include the parameters for hydrogen [2]. As discussed above, the role of the hydrogen atoms is to eliminate the effects of the unsaturated Si bonds on the surface of the wire, not to represent "real" hydrogen atoms. (In an experiment it is of course quite possible that hydrogen atoms, or other atoms, will attach themselves to these dangling bonds.)
There must naturally be an onsite element for hydrogen also:
Es_H = 0.999840*eV
onsite_H = SlaterKosterOnsiteParameters(element = Hydrogen,
angular_momenta = [0], # s-model for H
occupations = [1], # 1s
ionization_potential = [Es_H],
onsite_hartree_shift = [0]*eV
)
The scaling of the offsite parameters is completely analogous to those for Si:
h_si_sss = [ -3.999720*eV / (1.0+x)**2 for x in epsilon ] + [0.0*eV] h_si_ss1s = [ -1.697700*eV / (1.0+x)**2 for x in epsilon ] + [0.0*eV] h_si_sps = [ 4.251750*eV / (1.0+x)**2 for x in epsilon ] + [0.0*eV] h_si_sds = [ -2.105520*eV / (1.0+x)**2 for x in epsilon ] + [0.0*eV]
The only difference compared to before is that now there are two elements, which introduces keywords of the form "si_h_sss":
from NanoLanguage import *
'''
Generate a Slater-Koster basis set for silicon based on
T. B. Boykin et al., PRB 69, 115201 (2004)
http://dx.doi.org/10.1103/PhysRevB.69.115201
with hydrogen parameters, to remove effects of dangling bonds, from
Y. Zheng et al., IEEE Trans Electron Devices 52, 1097 (2005)
http://dx.doi.org/10.1109/TED.2005.848077
'''
from math import sqrt
# a is the cubic lattice constant
a = 5.430*Ang
# Nearest and second nearest neighbor distances
d_n1 = sqrt(3.0)/4.0*a
d_n2 = sqrt(2.0)/2.0*a
epsilon = numpy.linspace(-0.10, 0.10, 21)
distances = [ d_n1*(1.0+x) for x in epsilon ] + [ 0.5*(d_n1+d_n2) ]
si_si_sss = [ -1.95933*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_s1s1s = [ -4.24135*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_ss1s = [ -1.52230*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_sps = [ 3.02562*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_s1ps = [ 3.15565*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_sds = [ -2.28485*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_s1ds = [ -0.80993*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_pps = [ 4.10364*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_ppp = [ -1.51801*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_pds = [ -1.35554*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_pdp = [ 2.38479*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_dds = [ -1.68136*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_ddp = [ 2.58880*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
si_si_ddd = [ -1.81400*eV / (1.0+x)**2 for x in epsilon ] + [ 0.0*eV]
Es = -2.15168*eV
Ep = 4.22925*eV
Ed = 13.78950*eV
Es1 = 19.11650*eV
onsite = SlaterKosterOnsiteParameters(element = Silicon,
angular_momenta = [0,1,2,0], # sp3d5s*
occupations = [2,2,0,0], # 2s^2 2p^2
ionization_potential = [Es,Ep,Ed,Es1],
onsite_hartree_shift = [0,0,0,0]*eV
)
Es_H = 0.999840*eV
onsite_H = SlaterKosterOnsiteParameters(element = Hydrogen,
angular_momenta = [0], # s-model for H
occupations = [1], # 1s
ionization_potential = [Es_H],
onsite_hartree_shift = [0]*eV
)
h_si_sss = [ -3.999720*eV / (1.0+x)**2 for x in epsilon ] + [0.0*eV]
h_si_ss1s = [ -1.697700*eV / (1.0+x)**2 for x in epsilon ] + [0.0*eV]
h_si_sps = [ 4.251750*eV / (1.0+x)**2 for x in epsilon ] + [0.0*eV]
h_si_sds = [ -2.105520*eV / (1.0+x)**2 for x in epsilon ] + [0.0*eV]
basis_set = SlaterKosterTable(silicon = onsite,
hydrogen = onsite_H,
si_si_sss = zip(distances,si_si_sss),
si_si_s1s1s = zip(distances,si_si_s1s1s),
si_si_ss1s = zip(distances,si_si_ss1s),
si_si_sps = zip(distances,si_si_sps),
si_si_s1ps = zip(distances,si_si_s1ps),
si_si_sds = zip(distances,si_si_sds),
si_si_s1ds = zip(distances,si_si_s1ds),
si_si_pps = zip(distances,si_si_pps),
si_si_ppp = zip(distances,si_si_ppp),
si_si_pds = zip(distances,si_si_pds),
si_si_pdp = zip(distances,si_si_pdp),
si_si_dds = zip(distances,si_si_dds),
si_si_ddp = zip(distances,si_si_ddp),
si_si_ddd = zip(distances,si_si_ddd),
h_si_sss = zip(distances,h_si_sss),
h_si_ss1s = zip(distances,h_si_ss1s),
h_si_sps = zip(distances,h_si_sps),
h_si_sds = zip(distances,h_si_sds),
)
The script above for the Si basis set should be straightforward to modify for other elements. For instance, you can easily obtain a basis set for Ge using parameters from Ref. [3], or GaAs and all other III-V binary compounds by taking parameters from Ref. [4].
For precisely these two references the orbital labels etc are basically identical for all systems. If, on the other hand, you find a different model in an article, with quite different orbitals etc, there would be a relatively large risk to make a typo in the conversion process.
To this end, you may want to consider an alternative method to setup the Si basis set:
from NanoLanguage import *
element = Silicon
angular_momenta = [0,1,2,0] # sp3d5s*
occupations = [2,2,0,0] # 2s^2 2p^2
# a is the cubic lattice constant
a = 5.430*Ang
Es = -2.15168*eV
Ep = 4.22925*eV
Ed = 13.78950*eV
Es1 = 19.11650*eV
ionization_potentials = [Es,Ep,Ed,Es1]
offsite = {
"sss" : -1.95933*eV,
"s*s*s" : -4.24135*eV,
"ss*s" : -1.52230*eV,
"sps" : 3.02562*eV,
"s*ps" : 3.15565*eV,
"sds" : -2.28485*eV,
"s*ds" : -0.80993*eV,
"pps" : 4.10364*eV,
"ppp" : -1.51801*eV,
"pds" : -1.35554*eV,
"pdp" : 2.38479*eV,
"dds" : -1.68136*eV,
"ddp" : 2.58880*eV,
"ddd" : -1.81400*eV,
}
# ----------------------------------------------------------------------------
# Usually no need to change anything below
from math import sqrt
from string import lower
# Nearest neighbor distance
d_n1 = sqrt(3.0)/4.0*a
# Second-nearest neighbor distance
d_n2 = sqrt(2.0)/2.0*a
epsilon = numpy.linspace(-0.10, 0.10, 21)
distances = [ d_n1*(1.0+x) for x in epsilon ] + [ 0.5*(d_n1+d_n2) ]
# Element label (Si for Silicon, etc); use for the offsite keywords
e = lower(element.name())
# Generate a keyword argument for each offsite parameter in the dict
args = ""
for p in offsite.keys():
args += "%s_%s_%s = zip(distances, \
[ %f*eV /(1.0+x)**2 for x in epsilon ] + [ 0.0*eV]),\n" \
% (e,e,p.replace("*","1"),offsite[p])
onsite = SlaterKosterOnsiteParameters(element = element,
angular_momenta = angular_momenta,
occupations = occupations,
ionization_potential = ionization_potentials,
onsite_hartree_shift = [0,]*len(occupations)*eV
)
exec("""basis_set = SlaterKosterTable(%s = onsite,%s)""" % (e,args))
To some extent the beginning of the script is self-explanatory, and the parts which are below the dashed line are generic, and don't need to be modified, even if the model changes. The script relies on two primary techniques:
a dictionary to store the matrix elements, where the keys (string labels) correspond directly to the Hamiltonian matrix element label (here it is possible to use "*", it is later replaced by "1");
the use of the exec statement to execute code which has been assembled dynamically in the script itself - this demonstrates Python in its glory, as an interpretive language.
An obvious test case is to compute the band structure of graphene in a simple script.
# Silicon
lattice = FaceCenteredCubic(5.4306*Angstrom)
elements = [Silicon, Silicon]
fractional_coordinates = [[ 0. , 0. , 0. ],
[ 0.25 , 0.25 , 0.25 ]]
bulk_configuration = BulkConfiguration(
bravais_lattice=lattice,
elements=elements,
fractional_coordinates=fractional_coordinates
)
from SiBasisSet import basis_set
numerical_accuracy_parameters = NumericalAccuracyParameters(
k_point_sampling=(9, 9, 9),
)
calculator = SlaterKosterCalculator(
basis_set=basis_set,
numerical_accuracy_parameters=numerical_accuracy_parameters,
)
bulk_configuration.setCalculator(calculator)
bulk_configuration.update()
nlsave('silicon.nc', bulk_configuration)
# Bandstructure
bandstructure = Bandstructure(
configuration=bulk_configuration,
route=[['L', 'G', 'X', 'W', 'K', 'L', 'W', 'X', 'U'],['K', 'G']],
points_per_segment=100,
bands_above_fermi_level=All
)
nlsave('silicon.nc', bandstructure)
Download this script from the link just above, as well as
the script SiBasis.py
defining the Si basis set,
and save them in the same directory.
The script must be run from the command line, in this directory.
It finishes in a few seconds,
and produces a very accurate band structure of silicon.
|
|
Note |
|---|---|
|
The energy bands in Si are degenerate at the U and K points, and by the special way we construct the route in the script, the band structure plot makes a discontinuous break in the route between these points, i.e. it goes from X to U, then jumps to K and continues to Γ. |
A more interesting example, where you also can see another interesting feature in ATK is to compute the complex band structure of a silicon nanowire. Here you also need to utilize the hydrogen saturation of the surface, to avoid states in the band gap.
The first step is to build the configuration shown in the figure below. The process for this is described in detail in a separate tutorial on nanowires, which you can find at the Tutorial page of the QuantumWise website. Follow the instructions there to construct the narrow Si nanowire with a hydrogen terminated surface.
Send this structure to the Script Generator, and insert a New Calculator. Open the New Calculator and select the ATK-SE: Slater-Koster (device) calculator.
Select 9 k-points in the C direction. All other options can be left at their default values.
From Analysis, insert a ComplexBandstructure block. To cover the large band gap of such a narrow wire, set the energy range from -4 to 2 eV. It is also necessary to use a large number of energy points, to get a detailed plot of the complex band structure. Use 501 points.
Finally, set the output file name to something
suitable, like si_h_cbs.nc
and save the script as si_h_cbs.py.
Also download and save
the Si-H basis set script SiHBasis.py
in the same directory.
Open the saved script in a text editor; you may use the built-in one in VNL, or your own favorite code editor.
To use the basis set defined above, locate the line
basis_set = DFTBDirectory("cp2k/nonscc/")
and replace it by
from SiHBasisSet import basis_set
Also, remove the line
pair_potentials = DFTBDirectory("cp2k/nonscc/")
as well as the line referring to the pair_potentials variable
in the SlaterKosterCalculator setup:
calculator = SlaterKosterCalculator(
basis_set=basis_set,
pair_potentials=pair_potentials, ### remove this line
numerical_accuracy_parameters=numerical_accuracy_parameters,
)
The script you have now created needs to be run on the command
line, in the directory where SiHBasis.py resides.
It takes a few minutes to complete.
To inspect the results, locate the created file si_h_cbs.nc
in VNL and plot the complex band structure.
| |
|
|
| Pi model for carbon |
|
Bibliography |