Introduction¶
GRABNES (GRAphene-BN-Electronic Structure) allows you to calculate the electronic properties of several types of layered materials, including the obvious combinations of graphene and h-BN, hence the name that was given to this package.
One can build the systems for several basis systems using built-in methods which is useful as a first approximation when considering the rigid configurations. In reality, we often want to include lattice relaxation effects where we read in the structure given as output from LAMMPS (for example) as input for the GRABNES code. The go-to method in that case is ReadXYZ.
In this tutorial, we are going to first use the build-in methods to build the systam and calculate some observables for them while after we are going to focus on the ReadXYZ approach.
A few useful post-processing code-snippets are given at the end of the notebook.
[ ]:
Built-in methods¶
Graphene¶
Gendata_graphene.in allows you to build graphene systems. Run it using the grabnes executable in bin to make sure it runs correctly. - Modify it to build systems of different sizes, using CellSize. Check in the source code what this flag means. Get familiar with grep -r "SuperCell ."-style commands to find all the occurences of a certain flag inside the source code. - Modify the output files generate.pos or pos to your needs to visualize it using VESTA, VMD,
xcrysden or any other of your favorite tools.
Graphene with hBN, using effective model¶
Let’s switch our attention to graphene on hBN using the effective model introduced in [Jeil Jung, PRB 89, 205414 (2014)]. For this, we can use (some of) the following input flags: - MoirePotential .true. - MoireJeil .true. - MoireOffDiag .true. - MoireOnlyH0 .false. - MoireOnlyHZ .false. - MoireNoH0AndHZ .false. - MoireH0AndHZ .true. - MoirePotC0 0-0.01013 - MoirePotCz 0-0.00901 - MoirePotCab 0.01134 - MoirePotPhi0 1.510233401750693 -
MoirePotPhiz 0.147131255943122 - MoirePotPhiab 0.342084533390889 - Latticepercent 0-0.018181818181818
Tasks: - The value for the last flag works for a 55x55-sized system. Why is that? - Modify the value for a 56x56-sized system. - Why is there a negative sign in front of this value? - Plot the onsite energies for each of the atoms using a scatter plot where the color matches those energies. Confirm that the resulting moire pattern is periodic. - Calculate the DOS for a single moire unit cell using Lanczos recursion. Plot out the data in the resulting generate.DOS. Why does the plot
look spiky? - Run a convergence study on this DOS by increasing the SuperCell value. - Calculate the DOS using exact diagonalization. - Calculate the electronic bandstructure of this system. - Add a magnetic field to this system and repeat the previous steps. Why can’t you use exact diagonalization for this? Hint: read up on the Peierls phase substitution, see for instance Ref. [A. Cresti, Phys. Rev. B 103, 045402 (2021)].
Twisted bilayer graphene¶
A common way to describe a moire system is using the 4 indices defined by Hermann.
Tasks: - Create a twisted bilayer system using the 31 30 30 31 indices and visualize the structure. - Calculate the electronic bandstructure and visualize the flat band. - Calculate the DOS and optimize the supercell size to get clear signatures.
ReadXYZ¶
Twisted bilayer graphene¶
Repeat the twisted bilayer graphene calculations from the previous section using the relaxed structure given for download here.
Compare the results for the previous rigid configuration with the relaxed configuration observables.
Graphene on h-BN¶
Calculate the bandstructure for this system using the input file and structure given in the shared folder.
[ ]:
Code snippets¶
[4]:
import numpy as np
import matplotlib.pyplot as plt
import os
def latexify(fig_width=None, fig_height=None, columns=1):
"""Set up matplotlib's RC params for LaTeX plotting.
Call this before plotting a figure.
Parameters
----------
fig_width : float, optional, inches
fig_height : float, optional, inches
columns : {1, 2}
"""
# code adapted from http://www.scipy.org/Cookbook/Matplotlib/LaTeX_Examples
# Width and max height in inches for IEEE journals taken from
# computer.org/cms/Computer.org/Journal%20templates/transactions_art_guide.pdf
assert(columns in [1,2])
if fig_width is None:
fig_width = 3.39 if columns==1 else 6.9 # width in inches
if fig_height is None:
golden_mean = (sqrt(5)-1.0)/2.0 # Aesthetic ratio
fig_height = fig_width*golden_mean # height in inches
#MAX_HEIGHT_INCHES = 7.0
#if fig_height > MAX_HEIGHT_INCHES:
# print("WARNING: fig_height too large:" + fig_height +
# "so will reduce to" + MAX_HEIGHT_INCHES + "inches.")
# fig_height = MAX_HEIGHT_INCHES
params = {'backend': 'ps',
'text.latex.preamble': [r'\usepackage{gensymb}'],
'axes.labelsize': 8, # fontsize for x and y labels (was 10)
'axes.titlesize': 8,
'font.size': 8, # was 10
'legend.fontsize': 8, # was 10
'xtick.labelsize': 8,
'ytick.labelsize': 8,
'text.usetex': True,
'figure.figsize': [fig_width,fig_height],
'font.family': 'serif',
'figure.autolayout': True,
'lines.linewidth' : 1.0
}
matplotlib.rcParams.update(params)
def format_axes(ax):
SPINE_COLOR = 'gray'
for spine in ['top', 'right']:
ax.spines[spine].set_visible(False)
for spine in ['left', 'bottom']:
ax.spines[spine].set_color(SPINE_COLOR)
ax.spines[spine].set_linewidth(0.5)
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
for axis in [ax.xaxis, ax.yaxis]:
axis.set_tick_params(direction='out', color=SPINE_COLOR)
#q_sigma = c0 * 22.2 # [nm] hopping fitting parameter
#q_pi = a_cc * 22.2 # [nm] hopping fitting parameter
#r0 = 0.184 * a # [nm] hopping fitting parameter
from mpl_toolkits.axes_grid1 import make_axes_locatable
def getBandstructure(dir1, nPlotBands, EShift):
f = open(dir1+'/generate.bands.dat','r')
print(f)
i = 0
j = 0
l = 0
k = []
E = []
kVec = []
EVec = []
counterPlotBands = 0
for line in f:
i = i+1
a = line.split()
if i==10:
nbands = int(a[5])
print(nbands)
nk = int(a[7])
print(nbands,nk)
if i>16:
if len(a)==0:
j = j+0.5
#print(j)
else:
if (int(j)==int(nbands/2-nPlotBands/2 + counterPlotBands) and int(j) < int(nbands/2+nPlotBands/2)):
if l<nk:
k.append(float(a[0]))
E.append((float(a[1])-EShift))
l = l+1
if l==nk:
kVec.append(k)
EVec.append(E)
l = 0
counterPlotBands += 1
k = []
E = []
return kVec, EVec
def getSingleBand3(dir1, EShift, bandIndex):
f = open(dir1+'/generate.bands.dat','r')
i = 0
j = 0
l = 0
k = []
E = []
kVec = []
EVec = []
counterPlotBands = 0
for line in f:
i = i+1
a = line.split()
if i==10:
nbands = int(a[5])
nk = int(a[7])
#print(nbands,nk)
if i>14:
if len(a)==0:
j = j+0.5
#print(j)
else:
#print(nbands)
if (int(j)==int(nbands/2) + bandIndex):
if l<nk:
k.append(float(a[0]))
E.append(float(a[1])-EShift)
l = l+1
if l==nk:
kVec.append(k)
EVec.append(E)
l = 0
counterPlotBands += 1
k = []
E = []
return kVec, EVec
def getGap(dir1, EShift):
Emin = 10000000
Emax = -10000000
bandIndexVec = [-1]
for bandIndex in bandIndexVec:
kVec, EVec = getSingleBand3(dir1, EShift, bandIndex)
i = 0
if (np.max(EVec)>Emax):
Emax = np.max(EVec)
#if (np.min(EVec)<Emin):
# Emin= np.min(EVec)
bandIndexVec = [0]
for bandIndex in bandIndexVec:
kVec, EVec = getSingleBand3(dir1, EShift, bandIndex)
i = 0
if (np.min(EVec)<Emin):
Emin= np.min(EVec)
gap = (Emin - Emax)*1000
return gap
def getGap2(dir1, EShift):
Emin = 10000000
Emax = -10000000
bandIndexVec = [-3]
for bandIndex in bandIndexVec:
kVec, EVec = getSingleBand3(dir1, EShift, bandIndex)
i = 0
if (np.max(EVec)>Emax):
Emax = np.max(EVec)
#if (np.min(EVec)<Emin):
# Emin= np.min(EVec)
bandIndexVec = [-2]
for bandIndex in bandIndexVec:
kVec, EVec = getSingleBand3(dir1, EShift, bandIndex)
i = 0
if (np.min(EVec)<Emin):
Emin= np.min(EVec)
gap = (Emin - Emax)*1000
return gap
def getEnergyCutPlot(dir1, chosenEnergy, Nx, Ny, cutoffMin, cutoffMax, nEnergy):
os.chdir(dir1)
data = np.genfromtxt("generate.spectral",skip_header=4)
kx = data[:,0]
ky = data[:,1]
kz = data[:,2]
energy = (data[:,3]*3.5)-0.31
spectral = data[:,4]
cutoffMin = cutoff
cutoffMax = cutoff
i = 0
j = 0
eigVec = []
spectralVec = []
eig = np.zeros(nEnergy)
spect = np.zeros(nEnergy)
k = np.zeros(shape=(Nx*Ny,2))
for el1, el2, el3, el4 in zip(kx,ky,energy,spectral):
i = i+1
if (i==1):
k[j,0] = el1
k[j,1] = el2
if (i<=nEnergy):
eig[i-1] = el3
spect[i-1] = el4
if (i==nEnergy):
eigVec.append(eig)
spectralVec.append(spect)
i = 0
j = j+1
eig = np.zeros(nEnergy)
spect = np.zeros(nEnergy)
eigTrueVec = []
for i, kki in enumerate(k):
eigTrue = np.zeros(nEnergy)
eig = eigVec[i]
spec = spectralVec[i]
eigTrue[np.logical_and(np.logical_and(eig>chosenEnergy-cutoffMin,eig<chosenEnergy+cutoffMax),spec>0.0)] = 1.0
specSum = np.sum(spec[np.logical_and(np.logical_and(eig>chosenEnergy-cutoffMin,eig<chosenEnergy+cutoffMax),spec>0.0)])
if (eigTrue.any()==1):
eigTrueVec.append(specSum)
else:
eigTrueVec.append(0)
return k, eigTrueVec
K1 = np.array([1.70276024584813, 0.0000])
K1p = np.array([1.70243417150001, 3.332186294066186E-002])
def getEnergyCutPlot2(dir1, chosenEnergy, Nx, Ny, cutoffMin, cutoffMax, nEnergy):
os.chdir(dir1)
data = np.genfromtxt("generate.spectral")
kx = data[:,0]
ky = data[:,1]
kz = data[:,2]
energy = (data[:,3]*3.5)-0.31
spectral = data[:,4]
cutoffMin = cutoff
cutoffMax = cutoff
i = 0
j = 0
eigVec = []
spectralVec = []
eig = np.zeros(nEnergy)
spect = np.zeros(nEnergy)
k = np.zeros(shape=(Nx*Ny,2))
for el1, el2, el3, el4 in zip(kx,ky,energy,spectral):
i = i+1
if (i==1):
k[j,0] = el1
k[j,1] = el2
if (i<=nEnergy):
eig[i-1] = el3
spect[i-1] = el4
if (i==nEnergy):
eigVec.append(eig)
spectralVec.append(spect)
i = 0
j = j+1
eig = np.zeros(nEnergy)
spect = np.zeros(nEnergy)
eigTrueVec = []
for i, kki in enumerate(k):
eigTrue = np.zeros(nEnergy)
eig = eigVec[i]
spec = spectralVec[i]
eigTrue[np.logical_and(np.logical_and(eig>chosenEnergy-cutoffMin,eig<chosenEnergy+cutoffMax),spec>0.0)] = 1.0
specSum = np.sum(spec[np.logical_and(np.logical_and(eig>chosenEnergy-cutoffMin,eig<chosenEnergy+cutoffMax),spec>0.0)])
if (eigTrue.any()==1):
eigTrueVec.append(specSum)
else:
eigTrueVec.append(0)
return k, eigTrueVec
def getDOS(fileName):
data = np.genfromtxt(fileName)
xdata = data[:,0]
ydata = data[:,1]
return xdata, ydata
K1 = np.array([1.70276024584813, 0.0000])
K1p = np.array([1.70243417150001, 3.332186294066186E-002])
[ ]:
[ ]:
[ ]: