Twisted bilayer graphene (in-built methods)¶
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 Lancoz DOS and optimize the supercell size to get clear signatures. Note that the TB model we use includes a large number of neighbors, because of this, the memory requirements are much larger than in our previous session where we only include 3 neighbors in the GBN effective model calculations. For the sake
of this session, consider small cells (i.e. SuperCell 5 that can still be calculated easily on the mem128 machines). You can go up to SuperCell 20 if you use one of the large memory machines. - Compare the previous DOS with the DOS calculated by exact diagonalization. Reminder: the Lancos DOS has been renormalized by the number of atoms. To match it with the diagonalization DOS, one needs to multiply Lanczos DOS by the number of two times the number of atoms and divide it by the first
nearest neighbor hopping (check TB.Hopping value in your input file). The energy axis is also given in units of the first nearest neighbor hopping. - Add Anderson disorder to your system and see how it affects the DOS. - Add a magnetic field and compare the DOS in presence and absence of disorder. Can you explain what you see based on your knowledge of Landau levels?
tBG(read from external file, i.e. use ReadXYZ)¶
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.
Twisted bilayer graphene on hBN (bonus)¶
Try to combine the knowledge from previous session with the knowledge from current session to simulate the same observables for tBG on hBN. Note that for exact diagonalization calculations, a commensurate cell is required. You can avoid this complication when considering very large supercells for Lanczos where the periodic boundary effects become increasingly negligible.
[ ]:
Code snippets¶
Some routines¶
[1]:
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])
Plot the bandstructure¶
[ ]:
nPlotBands = 40
ncol = 1
nrow = 1
emin = -0.2
emax = 0.2
#emin = -0.002
#emax = 0.007
latexify(columns=1)#, fig_height=1.7)
g1 = gridspec.GridSpec(nrow, ncol)#, height_ratios=[6,1])
g1.update(wspace=0.15, hspace=0.0) # set the spacing between axes.
fig, ((ax0)) = plt.subplots(nrow, ncol, sharex='col', sharey='row')
ax0 = subplot(g1[0])
kpoint = 0.0
gpoint = 0.032230
mpoint = 0.060143
kppoint = 0.076259
factor = kppoint
max1_08 = 0.078838
max1_12 = 0.078838
EShift = 0.302
dir1 = "/Users/nihao/github/moire_course/grabnes2/outputFiles/tBG/"#/useReadXYZ/"
kVec, EVec = getBandstructure(dir1, nPlotBands, EShift)
i = 0
for el1, el2 in zip(kVec,EVec):
i = i+1
if (i==1):
#ax0.plot(el1, el2,color="C1",label=r"KC$_{VV10}$")
ax0.plot(el1, el2,color="C0",label=r"GAP$_{20}$")
else:
ax0.plot(el1, el2,color="C0")
Plot the onsite energies¶
[ ]:
os.chdir("/Users/nihao/github/moire_course/grabnes2/inputFiles/graphene/effectiveGBN/")
data = np.genfromtxt("generate.pos")
xdata = data[:,0][::2]
ydata = data[:,1][::2]
data = np.genfromtxt("generate.e")
zdata = data[:,0][::2]
plt.figure(figsize=(20,10))
plt.scatter(xdata, ydata, c=zdata,cmap="gist_ncar")
for DOS visualization:
os.chdir("/Users/nihao/github/moire_course/grabnes2/outputFiles/graphene/")
xdata, ydata = getDOS("generate.DOS")
plt.plot(xdata,ydata)
[ ]: