from envtb.vasp import eigenval
from envtb.vasp import procar
from envtb.wannier90 import w90hamiltonian
import numpy
from envtb.quantumcapacitance import electrostatics, quantumcapacitance
from matplotlib import pyplot
"""
This file contains functions for often used procedures.
They can also be considered as use cases.
It's also recommended to use them as starting point for your own
functions: type ?? after the function name and press enter,
then you can copy, directly modify and execute the code.
"""
def plot_vasp_bandstructure(eigenval_filename, plot_filename, output='save'):
"""
Plot the bandstructure contained in a VASP EIGENVAL file.
The output format is determined by the file ending of
plot_filename.
output: determines if the plot is written to a file ('save' - default
value) or displayed ('show')
"""
eigenvaldata = eigenval.EigenvalData(eigenval_filename)
kpoints = eigenvaldata.kpoints()
energies = eigenvaldata.energies()
plotter = w90hamiltonian.BandstructurePlot()
plotter.plot(kpoints, energies)
if output == 'save':
plotter.save(plot_filename)
if output == 'show':
plotter.show()
def plot_vasp_and_wannier90_bandstructure(eigenval_filename,
wannier90hr_filename,
poscar_filename,
plot_filename, usedorbitals='all',
usedhoppingcells_rings='all',
output='save'):
"""
Plot the bandstructure contained in a VASP EIGENVAL file and
a wannier90_hr.dat file.
usedorbitals: a list of used orbitals to use. Default is 'all'. Note:
this only makes sense if the selected orbitals don't interact with
other orbitals.
usedhoppingcells_rings: If you don't want to use all hopping parameters,
you can set the number of 'rings' surrounding the main cell here.
If it is a list (e.g. range(5)), several plots are created.
The default value is 'all'.
output: determines if the plot is written to a file ('save' -
default value) or displayed ('show')
"""
# TODO: "ring" is a stupid word
eigenvaldata = eigenval.EigenvalData(eigenval_filename)
vasp_kpoints = eigenvaldata.kpoints()
vasp_energies = eigenvaldata.energies()
w90ham = w90hamiltonian.Hamiltonian(wannier90hr_filename, poscar_filename)
# When usedhoppingcells_rings is not a list: create a list with 1 element
if not isinstance(usedhoppingcells_rings, list):
usedhoppingcells_rings = [usedhoppingcells_rings]
for ring in usedhoppingcells_rings:
if ring == 'all':
cells = 'all'
else:
cells = w90ham.unitcells_within_zone(ring, 'd', numpy.inf)
w90_energies = w90ham.bandstructure_data(vasp_kpoints, 'd',
usedhoppingcells=cells,
usedorbitals=usedorbitals)
plotter = w90hamiltonian.BandstructurePlot()
plotter.plot(vasp_kpoints, vasp_energies, 'r-')
plotter.plot(vasp_kpoints, w90_energies, 'b--')
if output == 'save':
# TODO: dateiname abschneiden, nicht vorstellen
plotter.save(str(ring)+"_"+plot_filename)
if output == 'show':
plotter.show()
def TopPzBandNrAtGamma(procar_filename, gnrwidth_rings, pzoffset=0):
"""
Reads a VASP PROCAR file from a zigzag GNR calculation and determines
the highest \pi* band at the Gamma point (the first point in the PROCAR
file, actually) that is relevant for the wannier90 calculation.
procar_filename: path to the PROCAR file.
gnrwidth_rings: number of benzene rings in transverse direction.
pzoffset (optional): if there are other pz-character bands at the
gamma point below the \pi* points, the function counts wrong. This
value is a manual correction for this problem: It assumes that
pzoffset more pz-like bands are below the highest
\pi*-band at the gamma point.
Return:
highestgoodpzband,energyatgammaofhighestgoodpzband
highestgoodpzband: Band number of the highest pz band at the Gamma point.
energyatgammaofhighestgoodpzband: Energy of this band at the Gamma point.
"""
procarData = procar.ProcarData(procar_filename)
chargedata = procarData.chargedata()
energydata = procarData.energydata()
# Nr of carbon atoms in a GNR of that width
nrofpzbands = gnrwidth_rings*2+2
# Charge data at Gamma point
gammapointdata = chargedata[0]
# Sum the pz charge for a particular band at the gamma point over all
# ions. Do that for all bands.
gammapointpzdata = \
[sum([ion[2] for ion in band]) for band in gammapointdata]
# Select band indices where there is pz charge at gamma
selectpzbands = \
[i for i in range(len(gammapointpzdata)) if gammapointpzdata[i] > 0.]
# Get band index of highest pz band at gamma point (index starting with 0,
# like always)
highestgoodpzband = selectpzbands[nrofpzbands-1+pzoffset]
# Energy of that band at gamma point
energyatgammaofhighestgoodpzband = energydata[0][highestgoodpzband]
return highestgoodpzband, energyatgammaofhighestgoodpzband
[docs]def plot_zigzag_graphene_nanoribbon_pz_bandstructure(
wannier90hr_graphene, poscarfile, wannier90woutfile, outcarfile,
width, output=None, usedhoppingcells_rings='all'):
"""
Plot the \pi bandstructure of a zigzag graphene nanoribbon based on a
wannier90 calculation of bulk graphene. The \pz orbitals have to be the
first two orbitals in the wannier90 file.
A .dat file with the bandstructure data is also saved as output.dat
(each line representing one k-point).
wannier90hr_graphene: path to the wannier90_hr.dat file containing the
graphene bulk calculation
poscarfile: path to the VASP POSCAR file of the graphene bulk calculation
outcarfile: path to the VASP OUTCAR file
width: width of the ribbon (number of rings).
output: path to the output image file. If None, nothing will be saved.
Default is None.
usedhoppingcells_rings: If you don't want to use all hopping parameters,
you can set the number of 'rings' surrounding the main cell here. If it
is a list (e.g. range(5)), several plots are created.
The default value is 'all'.
Return:
Hamiltonian: the generated GNR Hamiltonian.
data: bandstructure data.
path: path through reciprocal space.
"""
if width % 2 == 0:
unitcells = width/2+1
get_rid_of = 1
else:
unitcells = width/2+2
get_rid_of = 3
# When usedhoppingcells_rings is not a list: create a list with 1 element
if not isinstance(usedhoppingcells_rings, list):
usedhoppingcells_rings = [usedhoppingcells_rings]
ham = w90hamiltonian.Hamiltonian.from_file(
wannier90hr_graphene, poscarfile, wannier90woutfile, outcarfile)
for ring in usedhoppingcells_rings:
if ring == 'all':
cells = 'all'
else:
cells = ham.unitcells_within_zone(ring, 'd', numpy.inf)
ham2 = ham.create_supercell_hamiltonian(
[[0, 0, 0], [1, 0, 0]],
[[1, -1, 0], [1, 1, 0], [0, 0, 1]],
usedorbitals=(0, 1), usedhoppingcells=cells)
ham3 = ham2.create_supercell_hamiltonian(
[[0, i, 0] for i in range(unitcells)],
[[1, 0, 0], [0, unitcells, 0], [0, 0, 1]])
ham4 = ham3.create_modified_hamiltonian(
ham3.drop_dimension_from_cell_list(1),
usedorbitals=range(1, ham3.nrorbitals()-get_rid_of))
path = ham4.point_path([[0, 0, 0], [0.5, 0, 0]], 100)
if output is not None:
ham4.plot_bandstructure(path, str(ring)+"_"+output, 'd')
numpy.savetxt(str(ring)+"_"+output+'.dat', numpy.real(
data), fmt="%12.6G")
data = ham4.bandstructure_data(
path, 'd') # data is calculated twice, that's kind of stupid
return ham4, data, path
[docs]def plot_armchair_graphene_nanoribbon_pz_bandstructure(
wannier90hr_graphene, poscarfile, wannier90woutfile, width, output,
usedhoppingcells_rings='all'):
"""
Plot the \pi bandstructure of an armchair graphene nanoribbon based on a
wannier90 calculation of bulk graphene. The \pz orbitals have to be the
first two orbitals in the wannier90 file.
In this example function, width is the number of rings in transversal
direction. A .dat file with the bandstructure
data is also saved as output.dat (each line representing one k-point).
wannier90hr_graphene: path to the wannier90_hr.dat file containing the
graphene bulk calculation
poscarfile: path to the VASP POSCAR file of the graphene bulk calculation
width: width of the ribbon (number of rings). Must be an even number.
output: path to the output image file.
usedhoppingcells_rings: If you don't want to use all hopping parameters,
you can set the number of 'rings' surrounding the main cell here. If it
is a list (e.g. range(5)), several plots are created.
The default value is 'all'.
"""
unitcells = width
# When usedhoppingcells_rings is not a list: create a list with 1 element
if not isinstance(usedhoppingcells_rings, list):
usedhoppingcells_rings = [usedhoppingcells_rings]
ham = w90hamiltonian.Hamiltonian.from_file(
wannier90hr_graphene, poscarfile, wannier90woutfile)
for ring in usedhoppingcells_rings:
if ring == 'all':
cells = 'all'
else:
cells = ham.unitcells_within_zone(ring, 'd', numpy.inf)
ham2 = ham.create_supercell_hamiltonian(
[[0, 0, 0], [1, 0, 0]],
[[1, -1, 0], [1, 1, 0], [0, 0, 1]],
usedorbitals=(0, 1), usedhoppingcells=cells)
ham3 = ham2.create_supercell_hamiltonian(
[[i, 0, 0] for i in range(unitcells)],
[[unitcells, 0, 0], [0, 1, 0], [0, 0, 1]])
ham4 = ham3.create_modified_hamiltonian(
ham3.drop_dimension_from_cell_list(0))
path = ham4.point_path([[0, 0, 0], [0, 0.5, 0]], 100)
ham4.plot_bandstructure(path, str(ring)+"_"+output, 'd')
data = ham4.bandstructure_data(path, 'd')
numpy.savetxt(str(ring)+"_"+output+'.dat', numpy.real(
data), fmt="%12.6G")
[docs]def plot_zigzag_graphene_nanoribbon_pz_bandstructure_nn(
nnfile, width, output, length, magnetic_B=None):
"""
Plot the \pi bandstructure of a zigzag graphene nanoribbon based on a
n-th nearest neighbour parameterization of bulk graphene.
A .dat file with the bandstructure data is also saved as
output.dat (each line representing one k-point).
nnfile: path to the nearest-neighbour input file (see example files)
width: width of the ribbon (number of rings).
output: path to the output image file.
length: length in x-direction of the returned (not the bandstructure
plot!) unit cell (in rings)
Return: Hamiltonian of the unit cell with given length
"""
if width % 2 == 0:
unitcells = width/2+1
get_rid_of = 1
else:
unitcells = width/2+2
get_rid_of = 3
ham = w90hamiltonian.Hamiltonian.from_nth_nn_list(nnfile)
ham2 = ham.create_supercell_hamiltonian(
[[0, 0, 0], [1, 0, 0]],
[[1, -1, 0], [1, 1, 0], [0, 0, 1]])
ham3 = ham2.create_supercell_hamiltonian(
[[0, i, 0] for i in range(unitcells)],
[[1, 0, 0], [0, unitcells, 0], [0, 0, 1]])
ham4 = ham3.create_modified_hamiltonian(
ham3.drop_dimension_from_cell_list(1),
usedorbitals=range(1, ham3.nrorbitals()-get_rid_of),
magnetic_B=magnetic_B)
path = ham4.point_path([[-0.5, 0, 0], [0.5, 0, 0]], 100)
ham4.plot_bandstructure(path, output, 'd')
ham5 = ham4.create_supercell_hamiltonian([[i, 0, 0] for i in range(
length)], [[length, 0, 0], [0, 1, 0], [0, 0, 1]])
# data=ham4.bandstructure_data(path, 'd')
# numpy.savetxt(output+'.dat', numpy.real(data), fmt="%12.6G")
return ham5
[docs]def plot_armchair_graphene_nanoribbon_pz_bandstructure_nn(
nnfile, width, output, magnetic_B=None):
"""
Plot the \pi bandstructure of an armchair graphene nanoribbon based on
a n-th nearest neighbour parameterization of bulk graphene.
In this example function, width is the number of rings in transversal
direction. A .dat file with the bandstructure
data is also saved as output.dat (each line representing one k-point).
nnfile: path to the nearest-neighbour input file (see example files)
width: width of the ribbon (number of rings). Must be an even number.
output: path to the output image file.
"""
unitcells = width
ham = w90hamiltonian.Hamiltonian.from_nth_nn_list(nnfile)
ham2 = ham.create_supercell_hamiltonian(
[[0, 0, 0], [1, 0, 0]],
[[1, -1, 0], [1, 1, 0], [0, 0, 1]])
ham3 = ham2.create_supercell_hamiltonian(
[[i, 0, 0] for i in range(unitcells)],
[[unitcells, 0, 0], [0, 1, 0], [0, 0, 1]])
ham4 = ham3.create_modified_hamiltonian(
ham3.drop_dimension_from_cell_list(0),
magnetic_B=magnetic_B, gauge_B='landau_y')
path = ham4.point_path([[0, -0.5, 0], [0, 0.5, 0]], 100)
ham4.plot_bandstructure(path, output, 'd')
data = ham4.bandstructure_data(path, 'd')
numpy.savetxt(output+'.dat', numpy.real(data), fmt="%12.6G")
[docs]def PlotGraphenenthNNBandstructure(nnfile, outfile):
"""
Plot bandstructure of Graphene along the default hexagonal path.
nnfile: path to the nearest neighbour parameter file, e.g.
/path/to/envtb/exampledata/02_graphene_3rdnn/graphene3rdnnlist.dat
outfile: plot image output path, e.g. /tmp/bandstructure.png
"""
ham = w90hamiltonian.Hamiltonian.from_nth_nn_list(nnfile)
path = ham.standard_paths('hexagonal', 100)[2]
ham.plot_bandstructure(path, outfile, 'd')
[docs]def PlotGrapheneWannierBandstructure(
wannier90hr_graphene, poscarfile, wannier90woutfile, outfile):
"""
Plot bandstructure of Graphene along the default hexagonal path from
Wannier90 data.
You can use the files in
/path/to/envtb/exampledata/01_graphene_vasp_wannier90
wannier90hr_graphene: path to the wannier90_hr.dat file containing the
graphene bulk calculation
poscarfile: path to the VASP POSCAR file of the graphene bulk calculation
wannier90woutfile: path to the wannier90.wout file
output: path to the output image file.
"""
ham = w90hamiltonian.Hamiltonian.from_file(
wannier90hr_graphene, poscarfile, wannier90woutfile)
path = ham.standard_paths('hexagonal', 100)[2]
ham.plot_bandstructure(path, outfile, 'd')
[docs]def SimpleElectrostaticProblem():
"""
Calculates the electrostatic potential of a 1D stripe of 1nm width and
30 nm height with a gridsize of 1 nm. At 11 nm, the potential is fixed
to 8V.
Boundary conditions:
* left, right: Periodic boundary conditions (using PeriodicContainer).
* top: Neumann boundary condition, slope=1e18
* bottom: Dirichlet boundary condition, potential=0
(default boundary condition)
"""
breite = 1
hoehe = 30
gridsize = 1e-9
lapl = electrostatics.Laplacian2D2ndOrderWithMaterials(gridsize, gridsize)
rect = electrostatics.Rectangle(hoehe, breite, 1., lapl)
rect[hoehe-1, 0].neumannbc = (1e18, 'xb')
rect[10, 0].potential = 8
rect[0, 0].neumannbc = (0, 'xf')
cont = electrostatics.PeriodicContainer(rect, 'y')
solver, inhomogeneity = cont.lu_solver()
sol = solver(inhomogeneity)
pyplot.plot(sol)
print sol
[docs]def PotentialOfGluedRectangles2D():
"""
Two rectangles with capacitor plates, glued together with Container.
There are four capacitor plates. The boundary conditions
are the default BC.
The relative position of the two rectangles to each other is given by the
parameters position and offset of the connect() function (see
documentation).
"""
lapl = electrostatics.Laplacian2D2ndOrder(1, 1)
breite = 400
hoehe = 200
graphene_breite = 350
abstand = 20
rechteck = electrostatics.Rectangle(hoehe, breite, 1., lapl)
rechteck2 = electrostatics.Rectangle(300, 300, 1, lapl)
for x in range((breite-graphene_breite)/2, (breite+graphene_breite)/2):
rechteck[hoehe/2+5, x].potential = 5
rechteck[hoehe/2-5, x].potential = -5
for x in range(breite):
for y in range(hoehe/2-5, hoehe/2+5):
rechteck[y, x].epsilon = 1.
for x in range(100, 200):
rechteck2[x, 200].potential = 5
for x in range(150, 250):
rechteck2[x, 250].potential = 5
container = electrostatics.Container((rechteck, rechteck2))
container.connect(
rechteck, rechteck2, align='left', position='top', offset=(0, 150))
solver, inhomogeneity = container.lu_solver()
sol = solver(inhomogeneity)
pyplot.imshow(container.vector_to_datamatrix(sol)[0])
[docs]def PotentialOfSimpleConductor2D():
"""
Rectangle with 5V conducting plate at the bottom.
Boundary conditions:
* Left, right, top: default boundary condition 0V
* bottom: Dirichlet BC, 5V
"""
lapl = electrostatics.Laplacian2D2ndOrderWithMaterials(1e-9, 1e-9)
breite = 400
hoehe = 600
rechteck = electrostatics.Rectangle(hoehe, breite, 1., lapl)
for x in range(breite):
rechteck[hoehe-1, x].potential = 5
for x in range(breite):
for y in range(hoehe/2, hoehe):
rechteck[y, x].epsilon = 1.
container = electrostatics.Container((rechteck,))
solver, inhomogeneity = container.lu_solver()
sol = solver(inhomogeneity)
pyplot.imshow(container.vector_to_datamatrix(sol)[0])
[docs]def GrapheneQuantumCapacitance(equation='potential'):
"""
Calculate the quantum capacitance of Graphene on SiO2, including
temperature.
equation: 'charge' or 'potential', see documentation of
QuantumCapacitanceSolver.
Boundary conditions:
* left, right: periodic BC
* top: Neumann BC, slope=0
* bottom: backgate capacitor plate
Please read the code comments for further documentation.
Return:
voltages: The voltages where the capacitance is calculated (x values)
capacitance: the quantum capacitance (y(x) values)
"""
# Set system parameters
gridsize = 1e-9
hoehe = 600
breite = 1
backgatevoltage = 0
graphenepos = 299
temperature = 300
sio2 = 3.9
vstart = -60
vend = 60
dv = 0.5
# Finite difference operator
lapl = electrostatics.Laplacian2D2ndOrderWithMaterials(gridsize, gridsize)
# Create Rectangle
periodicrect = electrostatics.Rectangle(hoehe, breite, 1., lapl)
# Set boundary condition at the top
for y in range(breite):
periodicrect[0, y].neumannbc = (0, 'xf')
# Create list of backgate and GNR elements
backgateelements = [periodicrect[hoehe-1, y] for y in range(breite)]
grapheneelements = [periodicrect[graphenepos, y] for y in range(breite)]
# Set initial backgate element boundary condition (not necessary)
for element in backgateelements:
element.potential = backgatevoltage
# Create D(E,T) function.
Ef_dependence_function = quantumcapacitance.BulkGrapheneWithTemperature(
temperature, gridsize).Ef_interp
Ef_dependence_function_2 = quantumcapacitance.BulkGrapheneWithTemperature(
temperature, gridsize, interpolation=False).Q
# Set electrochemical potential and D(E,T) for the GNR elements
for element in grapheneelements:
element.potential = 0
element.fermi_energy = 0
element.fermi_energy_charge_dependence = Ef_dependence_function
element.charge_fermi_energy_dependence = Ef_dependence_function_2
# Set dielectric material
for x in range(graphenepos, hoehe):
for y in range(breite):
periodicrect[x, y].epsilon = sio2
# Create periodic container
percont = electrostatics.PeriodicContainer(periodicrect, 'y')
# Invert discretization matrix
solver, inhomogeneity = percont.lu_solver()
# Create QuantumCapacitanceSolver object.
# Depending on equation, either fermi_energy_charge_dependence or
# charge_fermi_energy_dependence will be used.
# You don't have to set the other one.
qcsolver = quantumcapacitance.QuantumCapacitanceSolver(
percont, solver, grapheneelements, lapl, equation=equation)
# Refresh basisvectors for calculation. Necessary once at the beginning.
qcsolver.refresh_basisvecs()
# Create volgate list
voltages = numpy.arange(vstart, vend, dv)
charges = []
# Loop over voltages
for v in voltages:
# print v
# Set backgate elements to voltage
for elem in backgateelements:
elem.potential = v
# Check change of environment
qcsolver.refresh_environment_contrib()
# Solve quantum capacitance problem & set potential property of
# elements
qcsolution = qcsolver.solve()
# Create new inhomogeneity because potential has changed
inhom = percont.createinhomogeneity()
# Solve system with new inhomogeneity
sol = solver(inhom)
# Save the charge configuration in the GNR
charges.append(percont.charge(sol, lapl, grapheneelements))
# Sum over charge and take derivative = capacitance
totalcharge = numpy.array([sum(x) for x in charges])
capacitance = (totalcharge[2:]-totalcharge[
:-2])/len(grapheneelements)*gridsize/(2*dv)
# Plot result
ax = pyplot.gca()
ax.set_title('Quantum capacitance of Graphene on SiO2')
ax.set_xlabel('Backgate voltage [V]')
ax.set_ylabel('GNR capacitance [$10^{-6} F/m^2$]')
pyplot.plot(voltages[1:-1], 1e6*capacitance)
return voltages[1:-1], capacitance