r"""
graphene_flake
~~~~~~~~~~~~~~
The **graphene_flake** package models a tight-binding graphene flake submitted to
inhomenehous strain leading to a uniform pseudo magnetic field and pseudo Landau
levels (pLLs).
**graphene_flake** is available at https://github.com/cpoli/graphene_flake
"""
import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as plt
from numpy import arange
import scipy.sparse as sparse
from math import sqrt
import os
[docs]class graphene_flake:
r"""
This is the mother class of the **graphene_flake** package which solve
the eigenvalue problem, plot the spectrum, and get the Hofstater butterfly.
The methods of this class are generic and works for any kind of flake.
Example usage::
# |
# B
# / \
# A A
#
# with the three hoppings:
# t_1 |
# t_2 /
# t_3 \
:param t_1: Hopping parameter t_1.
:param t_2: Hopping parameter t_2.
:param t_3: Hopping parameter t_3.
:param ext: Format of the figures to be saved.
"""
def __init__(self, t_1, t_2, t_3, ext):
if not isinstance(t_1, (int, float)) or t_1 <= 0.:
raise Exception('\n\nArgument t_1 not a positive number.\n')
if not isinstance(t_2, (int, float)) or t_2 <= 0.:
raise Exception('\n\nArgument t_2 not a positive number.\n')
if not isinstance(t_3, (int, float)) or t_3 == 0.:
raise Exception('\n\Argument t_3 not a positive number.\n')
self.t_1 = t_1
self.t_2 = t_2
self.t_3 = t_3
self.sites = -1
self.sites_a = -1
self.sites_b = -1
self.dx = sqrt(3.)/2.
self.dy = 0.5
self.strain = 0.
self.ind_max_x = 0
self.ind_max_y = 0
self.e_min = 0.
self.e_max = 0.
self.vortex = False
self.vacancy = False
self.vacancy_ind = -1
self.ext = ext
self.ind_a = np.array([])
self.ind_b = np.array([])
self.sites_x = np.array([])
self.sites_tot_x = np.array([])
self.bonds_x = np.array([])
self.bonds_tot_x = np.array([])
self.bonds_y = np.array([])
self.bonds_tot_y = np.array([])
self.coor_x = np.array([])
self.coor_y = np.array([])
self.coor_hop_x = np.array([])
self.coor_hop_y = np.array([])
self.hop_x = []
self.hop_y = []
self.eig_ener = np.array([])
self.eig_vec = np.array([])
self.eig_pola = np.array([])
self.energy_selec = np.array([])
self.density_selec = np.zeros([])
self.check_dir()
[docs] def get_indices(self):
"""
Get the indices associated with A and B sublattices.
"""
print('get_indices', self.sites_tot_x[-1])
self.ind_a = -np.ones(self.sites_tot_x[-1], int)
self.ind_a = -np.ones(self.sites_tot_x[-1], int)
self.ind_b = -np.ones(self.sites_tot_x[-1], int)
for i in range(self.rows_a):
self.ind_a[self.sites_tot_x[2*i]: self.sites_tot_x[2*i+1]] = \
range(self.sites_tot_x[2*i], self.sites_tot_x[2*i+1])
for i in range(self.rows_b):
self.ind_b[self.sites_tot_x[2*i+1]: self.sites_tot_x[2*i+2]] = \
range(self.sites_tot_x[2*i+1], self.sites_tot_x[2*i+2])
ind_b = np.flatnonzero(self.ind_a==-1)
ind_a = np.flatnonzero(self.ind_b==-1)
self.ind_b = np.delete(self.ind_b, ind_a)
self.ind_a = np.delete(self.ind_a, ind_b)
[docs] def get_eig_lin_strain(self, strain, eig_vec=False, vortex=False, vacancy=False, row_vacancy=0):
r"""
Build up the Hamiltonian :math:`H` and solve the eigenvalue problem
:math:`H|\psi_n\rangle=E_n|\psi_n\rangle`.
:param strain: Default value 0. Strength of the strain.
:param eig_vec: Default value False. Get and store the eigenvectors.
:param vortex: Default value False. Create a vortex at the middle of flake.
.. note::
The vortex is created by changing the sign of the hopppings
:math:`t_1\rightarrow -t_1`, for half of the bonds in the row
at the middle of the flake.
"""
if not isinstance(strain, (int, float)):
raise Exception('\n\nArgument strain not a number.\n')
if not isinstance(eig_vec, bool):
raise Exception('\n\nArgument eig_vec not a boolean.\n')
if not isinstance(vortex, bool):
raise Exception('\n\Argument vortex not a boolean.\n')
self.strain = strain
self.vortex = vortex
self.vacancy = vacancy
self.row_vacancy = row_vacancy
# hoppings t_1 (along the y direction)
ind_ham_y_a = np.zeros(self.bonds_tot_y[self.ind_max_y])
ind_ham_y_b = np.zeros(self.bonds_tot_y[self.ind_max_y])
self.hop_y = np.zeros(self.bonds_tot_y[self.ind_max_y])
for i in range(self.ind_max_y):
ym_1 = 0.5*(self.coor_y[self.sites_tot_x[2*i+1]]
+ self.coor_y[self.sites_tot_x[2*i+2]])
ind_ham_y_a[self.bonds_tot_y[i]: self.bonds_tot_y[i+1]] = \
arange(self.sites_tot_x[2*i+2], self.sites_tot_x[2*i+3])
ind_ham_y_b[self.bonds_tot_y[i]: self.bonds_tot_y[i+1]] = \
arange(self.sites_tot_x[2*i+1], self.sites_tot_x[2*i+2])
self.hop_y[self.bonds_tot_y[i]:self.bonds_tot_y[i+1]] = \
self.t_1*(1.+0.5*strain*ym_1)
# hoppings t_2 and t_3 (along the x direction)
len_x = len(self.sites_tot_x)//2-1
ind_ham_x_a = np.zeros(self.bonds_tot_x[-1])
ind_ham_x_b = np.zeros(self.bonds_tot_x[-1])
self.hop_x = np.zeros(self.bonds_tot_x[-1])
for i in range(self.ind_max_x):
# AB row starting with a B site
if (self.coor_x[self.sites_tot_x[2*i+1]]-self.coor_x[self.sites_tot_x[2*i]]) < 0.:
xm_2 = 0.5*(self.coor_x[self.sites_tot_x[2*i]:self.sites_tot_x[2*i+1]]
+ self.coor_x[1+self.sites_tot_x[2*i+1]:self.sites_tot_x[2*i+2]])
xm_3 = 0.5*(self.coor_x[self.sites_tot_x[2*i]:self.sites_tot_x[2*i+1]]
+ self.coor_x[self.sites_tot_x[2*i+1]:self.sites_tot_x[2*i+2]-1])
ym = 0.5*(self.coor_y[self.sites_tot_x[2*i]]+self.coor_y[self.sites_tot_x[2*i+1]])
sc_t2 = -self.dx*xm_2-0.5*ym
sc_t3 = self.dx*xm_3-0.5*ym
ind_ham_x_a[self.bonds_tot_x[2*i]:self.bonds_tot_x[2*i+1]] = \
arange(0+self.sites_tot_x[2*i], self.sites_tot_x[2*i+1])
ind_ham_x_b[self.bonds_tot_x[2*i]:self.bonds_tot_x[2*i+1]] = \
arange(1+self.sites_tot_x[2*i+1], self.sites_tot_x[2*i+2])
self.hop_x[self.bonds_tot_x[2*i]:self.bonds_tot_x[2*i+1]] = \
self.t_2*(1.+0.5*strain*sc_t2)
ind_ham_x_a[self.bonds_tot_x[2*i+1]:self.bonds_tot_x[2*i+2]] = \
arange(self.sites_tot_x[2*i], self.sites_tot_x[2*i+1])
ind_ham_x_b[self.bonds_tot_x[2*i+1]:self.bonds_tot_x[2*i+2]] = \
arange(self.sites_tot_x[2*i+1], self.sites_tot_x[2*i+2]-1)
self.hop_x[self.bonds_tot_x[2*i+1]:self.bonds_tot_x[2*i+2]] = \
self.t_3*(1.+0.5*strain*sc_t3)
# AB row starting with a A site
else:
xm_2 = 0.5*(self.coor_x[self.sites_tot_x[2*i]:self.sites_tot_x[2*i+1]-1]
+ self.coor_x[self.sites_tot_x[2*i+1]:self.sites_tot_x[2*i+2]])
xm_3 = 0.5*(self.coor_x[1+self.sites_tot_x[2*i]:self.sites_tot_x[2*i+1]]
+ self.coor_x[self.sites_tot_x[2*i+1]:self.sites_tot_x[2*i+2]])
ym = 0.5*(self.coor_y[self.sites_tot_x[2*i]]+self.coor_y[self.sites_tot_x[2*i+1]])
sc_t2 = -self.dx*xm_2-0.5*ym
sc_t3 = self.dx*xm_3-0.5*ym
ind_ham_x_a[self.bonds_tot_x[2*i]: self.bonds_tot_x[2*i+1]] = \
arange(1+self.sites_tot_x[2*i], self.sites_tot_x[2*i+1])
ind_ham_x_b[self.bonds_tot_x[2*i]: self.bonds_tot_x[2*i+1]] = \
arange(0+self.sites_tot_x[2*i+1], self.sites_tot_x[2*i+2])
self.hop_x[self.bonds_tot_x[2*i]: self.bonds_tot_x[2*i+1]] = \
self.t_3*(1.+0.5*strain*sc_t3)
ind_ham_x_a[self.bonds_tot_x[2*i+1]: self.bonds_tot_x[2*i+2]] = \
arange(self.sites_tot_x[2*i], self.sites_tot_x[2*i+1]-1)
ind_ham_x_b[self.bonds_tot_x[2*i+1]: self.bonds_tot_x[2*i+2]] = \
arange(self.sites_tot_x[2*i+1], self.sites_tot_x[2*i+2])
self.hop_x[self.bonds_tot_x[2*i+1]: self.bonds_tot_x[2*i+2]] = \
self.t_2*(1.+0.5*strain*sc_t2)
# hexagonal flake
'''
if vortex:
i_vortex = self.ind_max_y//2
self.hop_y[self.bonds_tot_y[i_vortex]:
self.bonds_tot_y[i_vortex]+self.bonds_y[i_vortex]//2+1] *= -1
'''
# triangular flake
if vortex:
i_vortex = self.ind_max_y//3
self.hop_y[self.bonds_tot_y[i_vortex]:
self.bonds_tot_y[i_vortex]+self.bonds_y[i_vortex]//2] *= -1
iA = np.concatenate([ind_ham_x_a, ind_ham_x_b, ind_ham_y_a, ind_ham_y_b])
iB = np.concatenate([ind_ham_x_b, ind_ham_x_a, ind_ham_y_b, ind_ham_y_a])
hop = np.concatenate([self.hop_x, self.hop_x, self.hop_y, self.hop_y])
ham = sparse.csr_matrix((hop, (iA, iB)), shape=(self.sites, self.sites))
ham = ham.toarray()
if vacancy:
ham = np.delete(ham, self.sites_tot_x[self.row_vacancy]+self.sites_x[self.row_vacancy+1]//2, 0)
ham = np.delete(ham, self.sites_tot_x[self.row_vacancy]+self.sites_x[self.row_vacancy+1]//2, 1)
mask = np.ones(self.sites_tot_x[-1], dtype=bool)
mask[self.sites_tot_x[self.row_vacancy]+self.sites_x[self.row_vacancy+1]//2] = False
self.coor_x = self.coor_x[mask]
self.coor_y = self.coor_y[mask]
self.sites_x[self.row_vacancy+1] -= 1
self.sites -= 1
for i in range(1, len(self.sites_x)):
self.sites_tot_x[i] = np.sum(self.sites_x[1:i+1])
self.get_indices()
print('solve the eigenvalue problem')
if eig_vec:
self.eig_ener, self.eig_vec = la.linalg.eigh(ham)
self.eig_pola = \
np.array([np.dot(self.eig_vec[self.ind_a, i], self.eig_vec[self.ind_a, i]) for i in range(self.sites)])
else:
self.eig_ener = la.linalg.eigvalsh(ham)
[docs] def get_states_selec(self, e_min, e_max):
r"""
Get, if any, the wave function support of the states between *e_min* and *e_max*.
"""
if not isinstance(e_min, float) or not isinstance(e_max, float):
raise Exception('\n\nArguments e_min and e_max not float.\n')
if not self.eig_vec.any():
raise Exception('\n\nRun the method get_eig_lin_strain() with eig_vec=True.\n')
self.e_min, self.e_max = e_min, e_max
self.density_selec = np.zeros(len(self.eig_ener), float)
index = np.array(((self.eig_ener > e_min) & (self.eig_ener < e_max)).nonzero())
index = np.argwhere((self.eig_ener>e_min) & (self.eig_ener<e_max))
index = np.ravel(index)
if index.any():
m = len(index)
print('GET {} STATE(S) between {} and {}:'.format(m, e_min, e_max))
self.energy_selec = np.empty(m)
for i in range(m):
self.density_selec += np.abs(self.eig_vec[:, index[i]])**2
self.density_selec /= m
print('A sublattice polarization ', np.sum(self.density_selec[self.ind_a]))
print('B sublattice polarization ', np.sum(self.density_selec[self.ind_b]))
else:
print('NO STATES in the selected range.')
[docs] def get_butterfly(self, nbr_points=100, ener_lim=[-2., 2.], vortex=False, save=True):
strain_max_neg, strain_max_pos = self.get_lin_strain_lims()
strain_vec = np.linspace(strain_max_neg, strain_max_pos, nbr_points)
self.get_eig_lin_strain(strain=0., eig_vec=False, vortex=vortex)
ind_en = np.argwhere((self.eig_ener>ener_lim[0]) & (self.eig_ener<ener_lim[1]))
ind_en = np.ravel(ind_en)
mat_ener = np.empty((len(ind_en), nbr_points))
for i, s in enumerate(strain_vec):
print(i, s)
self.get_eig_lin_strain(strain=s, eig_vec=False, vortex=vortex)
mat_ener[:, i] = self.eig_ener[ind_en]
fig, ax = plt.subplots()
fig.canvas.set_window_title('BUTTERFLY')
ax.set_xlabel(r'$\beta/\beta_{max}$', fontsize=20)
ax.set_ylabel('$E$', fontsize=20)
ax.set_xticks((strain_max_pos, 0.5*strain_max_pos, 0.,
-0.5*strain_max_pos, -strain_max_pos))
ax.set_xticklabels(('-1', '-1/2', '0', '1/2', '1'))
ax.set_xlim([strain_vec[0], strain_vec[-1]])
ax.set_ylim(ener_lim)
for i in range(len(ind_en)):
plt.plot(strain_vec, mat_ener[i, :], color='#00008B')
fig.set_tight_layout(True)
plt.draw()
if save:
name_fig = 'N'+str(self.sites)
if self.vortex:
name_fig += 'vor'
fig.savefig('figs/butterfly_'+name_fig+'.png', format='png')
[docs] def plt_spec(self, nbr_bins=101, ener_lim=[-1., 1.], ms=10, plt_pola=False, save=False):
"""
Plot the spectrum in two different ways:
* energy level staircase.
* number of states (histogram).
:param nbr_bins: Default value 101. Number of bins of the histogram.
:param ener_lim: Default value [-1.5, 1.5]. list of the energy min and max.
:param ms: Default value 10. Size of the markers.
:param save: Default value False. Save the two figures.
"""
if not isinstance(nbr_bins, int) or nbr_bins<=0:
raise Exception('\n\nArgument nbr_bins not a positive integers.\n')
if not isinstance(ener_lim, list) or len(ener_lim) != 2:
raise Exception('\n\nArgument ener_lim not a list of dim 2.\n')
if not isinstance(save, bool):
raise Exception('\n\Argument save not a boolean.\n')
fig, ax1 = plt.subplots()
fig.canvas.set_window_title('SPECTRUM')
ax1.set_xlabel('$j/N$', fontsize=20)
ax1.set_ylabel('E/t', fontsize=20)
ind_min = np.argmin(np.abs(self.eig_ener-ener_lim[0]))
ind_max = np.argmin(np.abs(self.eig_ener-ener_lim[1]))
ind_en = np.argwhere((self.eig_ener>ener_lim[0]) & (self.eig_ener<ener_lim[1]))
ind_en = np.ravel(ind_en)
en = self.eig_ener[ind_en]
ind = np.arange(self.sites)-self.sites//2
ind = ind[ind_en]
ax1.set_ylim([en[0], en[-1]])
ax1.set_xlim([ind[0], ind[-1]])
ax1.plot(ind, en, '.', color='#00008B', ms=ms)
if self.eig_pola.any() and plt_pola:
for label in ax1.get_yticklabels():
label.set_color('blue')
ax1.set_ylabel('E/t', fontsize=20, color='blue')
ax2 = ax1.twinx()
ax2.plot(ind, self.eig_pola[ind_en], '.r', ms=6)
ax2.set_yticks([0, 0.5, 1])
ax2.set_ylabel(r'$\langle A|A\rangle$', fontsize=20, color='red')
ax2.set_ylim([-0.1, 1.1])
ax2.set_xlim([ind[0], ind[-1]])
for label in ax2.get_yticklabels():
label.set_color('red')
if save:
name_fig = 'N'+str(self.sites)+'_'+str(self.strain).replace('.', ',')
if self.vortex:
name_fig += 'vor'
fig.savefig('figs/spec_'+name_fig+'.png', format='png')
fig, ax = plt.subplots()
fig.canvas.set_window_title('SPECTRUM')
n, bins, patches = plt.hist(en, bins=nbr_bins, color='#00008B')
ax.set_xlabel('$E$', fontsize=20)
ax.set_ylabel('number of states', fontsize=20)
ax.set_ylim([0, np.max(n)])
ax.set_xlim(ener_lim)
if save:
name_fig = 'N'+str(self.sites)+'_'+str(self.strain).replace('.', ',')
if self.vortex:
name_fig += 'vor'
fig.savefig('figs/spec_hist'+name_fig+'.png', format='png')
plt.draw()
[docs] def plt_states_selec(self, lw=1, coef=100000., fontsize=20,
plt_sites=False, plt_bonds=False, plt_axis=False, save=False):
"""
Plot the sum of the probability densities satisfying the condition
imposed by the method *get_states_cond*.
:param lw: Default value 1. Line width of the bonds.
:param coef: Default value 100000. Increase the area of the circles by a factor *coef*.
:param fontsize: Default value 20. Size of the font used for the axis labels and the title.
:param plt_bonds: Default value False. Plot the bonds.
:param plt_bonds: Default value False. Plot the sites.
:param save: Default value False. Save the figures.
"""
if not isinstance(coef, (int, float)) or coef<=0:
print(isinstance(coef, float), coef<=0)
raise Exception('\n\nArgument coef not a positive float.\n')
if not isinstance(fontsize, int) or fontsize<=0:
raise Exception('\n\nArgument fontsize not a positive integers.\n')
if not isinstance(plt_sites, bool):
raise Exception('\n\Argument plt_sites not a boolean.\n')
if not isinstance(plt_bonds, bool):
raise Exception('\n\Argument plt_bonds not a boolean.\n')
if not isinstance(save, bool):
raise Exception('\n\Argument save not a boolean.\n')
if self.density_selec.any():
fig, ax = plt.subplots()
fig.canvas.set_window_title('Sublattice polarized state(s)')
ax.set_xlabel('$x$', fontsize=fontsize)
ax.set_ylabel('$y$', fontsize=fontsize)
plt.scatter(self.coor_x[self.ind_a], self.coor_y[self.ind_a], c='b',
s=coef*self.density_selec[self.ind_a], alpha=0.5)
plt.scatter(self.coor_x[self.ind_b], self.coor_y[self.ind_b], c='r',
s=coef*self.density_selec[self.ind_b], alpha=0.5)
ax.set_aspect('equal')
ax.set_ylim([np.min(self.coor_y)-4., np.max(self.coor_y)+4.])
if not plt_axis:
ax.axis('off')
if plt_sites:
plt.scatter(self.coor_x, self.coor_y, s=5, c='k')
if plt_bonds:
self.plt_bonds(self.coor_x, self.coor_y, lw=lw)
plt.draw()
if save:
name_fig = 'N'+str(self.sites)+'_'+str(self.strain).replace('.', ',')
name_fig += '_'+str(self.e_min).replace('.', ',')+'_'+str(self.e_max).replace('.', ',')
if self.vortex:
name_fig += 'vor'
fig.savefig('figs/states_'+name_fig+'.'+self.ext, format=self.ext)
else:
print('NO SELECTED STATES\n')
[docs] def plt_lattice(self, lw=1, s=100, plt_sites=True, plt_bonds=True, plt_axis=False, save=False):
"""
Plot the lattice. Blue points correspond to the A sublattice.
Red points correspond the the B sublattice.
:param lw: Default value 1. Line width of the bonds.
:param s: Default value 100. Area of each marker.
:param plt_sites: Default True. Plot the sites.
:param plt_bonds: Default True. Plot the bonds.
:param plt_axis: Default False. Plot the axis.
:param save: Default False. Save the figure.
"""
if not isinstance(s, int) or s<=0:
raise Exception('\n\nArgument s not a positive integer.\n')
if not isinstance(plt_sites, bool):
raise Exception('\n\Argument plt_sites not a boolean.\n')
if not isinstance(plt_bonds, bool):
raise Exception('\n\Argument plt_bonds not a boolean.\n')
if not isinstance(plt_axis, bool):
raise Exception('\n\Argument plt_axis not a boolean.\n')
if not isinstance(save, bool):
raise Exception('\n\Argument save not a boolean.\n')
fig, ax = plt.subplots()
fig.canvas.set_window_title('LATTICE')
ax.set_xlim([np.amin(self.coor_x)-2, np.amax(self.coor_x)+2])
ax.set_ylim([self.coor_y[0]-2, np.amax(self.coor_y)+2])
ax.set_aspect('equal')
if not plt_axis:
ax.axis('off')
if plt_bonds:
self.plt_bonds(self.coor_x, self.coor_y, lw=lw)
print('self.vacancy', self.vacancy)
if plt_sites:
plt.scatter(self.coor_x[self.ind_a], self.coor_y[self.ind_a], s=s, c='b')
plt.scatter(self.coor_x[self.ind_b], self.coor_y[self.ind_b], s=s, c='r')
if save:
name_fig = 'N'+str(self.sites)+'_'+str(self.strain).replace('.', ',')
fig.savefig('figs/lattice_'+name_fig+'.png', format='png')
[docs] def plt_lattice_hop(self, lw=1, s=100, real_mw=False, plt_axis=False,
plt_sites=True, plt_bonds=True, save=True):
r"""
Plot the lattice in hopping space.
Blue points correspond to the A sublattice.
Red points correspond the the B sublattice.
:param lw: Default value 1. Line width of the bonds.
:param s: Default value 100. Area of each marker.
:param real_mw: Default value False. Plot the lattice in real space considering
the microwave regime (:math:`t\approx 10` and :math:`d\approx 10mm`)
:param plt_sites: Default False. Plot the sites.
:param plt_bonds: Default False. Plot the bonds.
:param plt_axis: Default False. Plot the axis.
:param save: Default False. Save the figure.
"""
if not isinstance(s, int) or s<=0:
raise Exception('\n\nArgument s not a positive integer.\n')
if not isinstance(real_mw, bool):
raise Exception('\n\Argument real_mw not a boolean.\n')
if not isinstance(plt_sites, bool):
raise Exception('\n\Argument plt_sites not a boolean.\n')
if not isinstance(plt_bonds, bool):
raise Exception('\n\Argument plt_bonds not a boolean.\n')
if not isinstance(plt_axis, bool):
raise Exception('\n\Argument plt_axis not a boolean.\n')
if not isinstance(save, bool):
raise Exception('\n\Argument save not a boolean.\n')
if real_mw:
hop_plt_x = hop_to_real_mw(np.abs(self.hop_x))
hop_plt_y = hop_to_real_mw(np.abs(self.hop_y))
else:
hop_plt_x = np.abs(self.hop_x)
hop_plt_y = np.abs(self.hop_y)
self.get_coor_lin_strain(hop_plt_x, hop_plt_y)
fig, ax = plt.subplots()
ax.set_aspect('equal')
if not plt_axis:
ax.axis('off')
if real_mw:
fig.canvas.set_window_title('REAL SPACE')
else:
fig.canvas.set_window_title('HOPPINGS SPACE')
ax.set_xlim([np.amin(self.coor_hop_x)-2., np.amax(self.coor_hop_x)+2.])
ax.set_ylim([np.amin(self.coor_hop_y)-2., np.amax(self.coor_hop_y)+2.])
if plt_sites:
plt.scatter(self.coor_hop_x[self.ind_a], self.coor_hop_y[self.ind_a], s=s, c='b')
plt.scatter(self.coor_hop_x[self.ind_b], self.coor_hop_y[self.ind_b], s=s, c='r')
if plt_bonds:
self.plt_bonds(self.coor_hop_x, self.coor_hop_y, lw=lw)
if save:
file_name = 'figs/coor_N{}_n{}_m{}_beta' \
.format(self.sites, self.sites_base, self.sites_cut) \
+ str(self.strain).replace('.', ',') + '.txt'
np.savetxt(file_name, np.c_[self.coor_x, self.coor_y])
name_fig = 'N'+str(self.sites)+'_'+str(self.strain).replace('.', ',')
fig.savefig('figs/lattice_hop_'+name_fig+'.png', format='png')
plt.draw()
[docs] def plt_bonds(self, coor_x, coor_y, lw=1):
"""
Plot the bonds of the lattice.
"""
# plot bonds along x
for i in range(0, self.size_array_x-2, 2):
if coor_x[self.sites_tot_x[i+1]] > coor_x[self.sites_tot_x[i]]:
plt.plot([coor_x[self.sites_tot_x[i]: self.sites_tot_x[i+1]-1],
coor_x[self.sites_tot_x[i+1]: self.sites_tot_x[i+2]]],
[coor_y[self.sites_tot_x[i]: self.sites_tot_x[i+1]-1],
coor_y[self.sites_tot_x[i+1] :self.sites_tot_x[i+2]]],
'k', lw=lw) # t2
plt.plot([coor_x[self.sites_tot_x[i]+1: self.sites_tot_x[i+1]],
coor_x[self.sites_tot_x[i+1]: self.sites_tot_x[i+2]]],
[coor_y[self.sites_tot_x[i]+1: self.sites_tot_x[i+1]],
coor_y[self.sites_tot_x[i+1]: self.sites_tot_x[i+2]]],
'k', lw=lw) # t3
else:
plt.plot([coor_x[self.sites_tot_x[i]: self.sites_tot_x[i+1]],
coor_x[self.sites_tot_x[i+1]+1: self.sites_tot_x[i+2]]],
[coor_y[self.sites_tot_x[i]: self.sites_tot_x[i+1]],
coor_y[self.sites_tot_x[i+1]+1: self.sites_tot_x[i+2]]],
'k', lw=lw) # t2
plt.plot([coor_x[self.sites_tot_x[i]: self.sites_tot_x[i+1]],
coor_x[self.sites_tot_x[i+1]: self.sites_tot_x[i+2]-1]],
[coor_y[self.sites_tot_x[i]: self.sites_tot_x[i+1]],
coor_y[self.sites_tot_x[i+1]: self.sites_tot_x[i+2]-1]],
'k', lw=lw) # t3
# plot bonds along y
for i in range(1, self.size_array_y, 2):
plt.plot([coor_x[self.sites_tot_x[i]: self.sites_tot_x[i+1]],
coor_x[self.sites_tot_x[i+1]: self.sites_tot_x[i+2]]],
[coor_y[self.sites_tot_x[i]: self.sites_tot_x[i+1]],
coor_y[self.sites_tot_x[i+1]: self.sites_tot_x[i+2]]],
'k', lw=lw) # t1
[docs] def check_dir(self):
"""
Create if the directory to store the figures exists.
"""
if not os.path.exists('figs'):
os.makedirs('figs')
[docs] def plt_show(self):
"""
Duplicate the matplotlib method *show* to avoid importing
matplotlib to plot on screen the figures generated
by the **graphene_flake** class.
"""
plt.show()
[docs]class tri_to_hexa_zz(graphene_flake):
r"""
This is a child class of **graphene_flake**. **tri_to_hexa_zz** models
a triangular flake and the transition from triangular to hexagonal flake.
The transition is obtained but cutting the edges of the triangular flake.
Example usage::
# Triangular flake with one 3 sites at the base:
# A 8
# | |
# B 7
# / \ / \
# A A 5 6
# | | | |
# B B 3 4
# / \ / \ / \ / \
# A A A 0 1 2
#
#
# Triangular flake with one 3 sites at the base and 1 cut:
# B 7
# / \ / \
# A A 5 6
# | | | |
# B B 3 4
# \ / \ /
# A 1
#
:param sites_base: Number of A sites at the base of the triangle. Min value 3.
:param sites_cut: Default value 0. Number of cuts.
* triangle: sites_cut = 0.
* hexagon 0 < sites_cut :math:`\le` (sites_base+2)/3.
* regular hexagon: sites_cut = (sites_base+2)/3.
.. note::
3*sites_cut*-2 :math:`\le` *sites_base*.
:param t_1: Default value 1. Hopping parameter t_1.
:param t_2: Default value 1. Hopping parameter t_2.
:param t_3: Default value 1. Hopping parameter t_2.
:param ext: Default value 'png'. Format of the saved figures.
"""
def __init__(self, sites_base, sites_cut=0, t_1=1., t_2=1., t_3=1., ext='png'):
graphene_flake.__init__(self, t_1, t_2, t_3, ext=ext)
if not isinstance(sites_base, int):
raise Exception('\n\nArgument sites_base not a positive integers >= 3.\n')
if not isinstance(sites_cut, int) or sites_cut<0 or sites_cut>((sites_base+2)//3):
raise Exception('\n\nArgument sites_cut not a positive integers'
'between 0 and (sites_base+2)/3.\n')
sites_dif = sites_base-sites_cut
self.sites_base = sites_base
self.sites_cut = sites_cut
len_x = 2*sites_dif+5
self.bonds_x = np.zeros(len_x, int)
self.sites_x = np.zeros(len_x, int)
self.bonds_tot_x = np.zeros(len_x, int)
self.sites_tot_x = np.zeros(len_x, int)
self.bonds_x[1: 2*sites_cut+1: 2] = \
arange(sites_base-2*(sites_cut-1), sites_dif+2)
self.bonds_x[2: 2*sites_cut+1: 2] = \
arange(sites_base-2*(sites_cut-1), sites_dif+2)
self.bonds_x[2*sites_cut+1: len_x: 2] = \
arange(sites_dif+1, sites_cut-1, -1)
self.bonds_x[2*sites_cut+2: len_x : 2] = \
arange(sites_dif+1, sites_cut-1, -1)
if sites_cut == 0:
self.sites_x[1] = sites_base+2
self.sites_x[2] = sites_base+1
else:
self.sites_x[1] = sites_base
self.sites_x[2] = sites_base+1
self.sites_x[3: len_x: 2] = arange(sites_base+1, sites_cut, -1)
self.sites_x[4: len_x: 2] = arange(sites_base, sites_cut-1, -1)
# vacancy, if needed, the edge sites at the bottom of the flake
for i in range(2, sites_cut+1):
self.sites_x[1: 2*i] -= 2
self.bonds_y = np.zeros(sites_dif+2, int)
self.bonds_tot_y = np.zeros(sites_dif+2, int)
if sites_cut == 0:
self.bonds_y[1: sites_dif+2] = arange(sites_base+1, sites_cut, -1)
else:
self.bonds_y[1: 2+sites_cut-1] = \
arange(sites_base+1-2*(sites_cut-1), sites_base-sites_cut+3)
self.bonds_y[1+sites_cut-1:sites_dif+2] = \
arange(sites_base+1-(sites_cut-1), sites_cut, -1)
for i in range(1, sites_dif+2):
self.bonds_tot_y[i] = np.sum(self.bonds_y[1:i+1])
for i in range(1, len_x):
self.sites_tot_x[i] = np.sum(self.sites_x[1:i+1])
self.bonds_tot_x[i] = np.sum(self.bonds_x[1:i+1])
self.sites = self.sites_tot_x[len_x-1]
self.sites_a = sum(self.sites_x[0: 2*sites_dif+4: 2])
self.sites_b = sum(self.sites_x[1: 2*sites_dif+4: 2])
# get the coordinates
self.coor_x = np.zeros(self.sites)
self.coor_y = np.zeros(self.sites)
y = 0
for i in range(len_x-1):
self.coor_x[self.sites_tot_x[i]: self.sites_tot_x[i+1]] = \
self.dx*arange(1-self.sites_x[i+1], self.sites_x[i+1]+1./2, 2.)
self.coor_y[self.sites_tot_x[i]: self.sites_tot_x[i+1]] = y
y += .5
if i%2 == 1:
y += .5
center_y = sum(self.coor_y)/len(self.coor_y)
self.coor_y -= center_y
if sites_cut == 0:
self.ind_max_x = sites_base + 1
self.size_array_x = 2*sites_base + 4
else:
self.ind_max_x = sites_dif + 2
self.size_array_x = 2*sites_dif + 5
self.ind_max_y = sites_dif+1
self.size_array_y = 2*sites_dif+2
self.rows_a = len(self.sites_tot_x[0::2])-1
self.rows_b = len(self.sites_tot_x[1::2])
self.get_indices()
print('number of sites', self.sites_tot_x[-1])
[docs] def get_lin_strain_lims(self):
"""
Get the maximal negative and positive strain keeping the hoppings positive.
"""
y_m = 0.5*(self.coor_y[self.sites_tot_x[1]]+self.coor_y[self.sites_tot_x[2]])
strain_max_pos = -2./y_m - 1e-4
print(self.coor_y[self.sites_tot_x[1]], self.coor_y[self.sites_tot_x[2]])
sites_dif = self.sites_base - self.sites_cut
i = sites_dif+1
y_m = 0.5*(self.coor_y[self.sites_tot_x[2*i-1]]+self.coor_y[self.sites_tot_x[2*i]])
print(self.coor_y[self.sites_tot_x[2*i-1]], self.coor_y[self.sites_tot_x[2*i]])
strain_max_neg = -2./y_m + 1e-4
print('\nMaximal strain values keeping positive hoppings are between:')
print(' ', strain_max_neg, ' ', strain_max_pos, '\n')
return strain_max_neg, strain_max_pos
[docs] def get_coor_lin_strain(self, hop_plt_x, hop_plt_y):
"""
Get the lattice site positions in hopping space.
Method essentially used to be plot the lattice in hopping space.
:param hop_plt_x: array containing the hoppings (or distances) along x.
:param hop_plt_y: array containing the hoppings (or distances) along y.
"""
sites_dif = self.sites_base - self.sites_cut
hop_ord_x = np.zeros(self.bonds_tot_x[-1])
for i in range(self.ind_max_y+1):
hop_ord_x[1+self.bonds_tot_x[2*i]: self.bonds_tot_x[2*(i+1)]: 2] = \
hop_plt_x[self.bonds_tot_x[2*i]: self.bonds_tot_x[2*i+1]]
hop_ord_x[0+self.bonds_tot_x[2*i]: self.bonds_tot_x[2*(i+1)]: 2] = \
hop_plt_x[self.bonds_tot_x[2*i+1]: self.bonds_tot_x[2*i+2]]
coor_ord_x = np.zeros(self.sites)
coor_ord_y = np.zeros(self.sites)
k= 0
# if the number of AB sites along x decreases
l = 1
for j in range(self.sites_cut):
for i in range(1+self.sites_tot_x[2*j], self.sites_tot_x[2*j+2]):
coor_ord_x[i] = coor_ord_x[i-1] + self.dx*hop_ord_x[k]
coor_ord_y[i] = coor_ord_y[i-1] + (-1)**l*self.dy*hop_ord_x[k]
k += 1
l += 1
coor_ord_x[self.sites_tot_x[2*j+2]] = \
coor_ord_x[self.sites_tot_x[2*j]] - self.dx*hop_ord_x[k]
coor_ord_y[self.sites_tot_x[2*j+2]] = \
coor_ord_y[self.sites_tot_x[2*j]] + self.dy*hop_ord_x[k] \
+ hop_plt_y[self.bonds_tot_y[j]]
# the interface, if any.
if self.sites_cut > 0:
coor_ord_x[self.sites_tot_x[2*j+2]] = coor_ord_x[self.sites_tot_x[2*j]]
coor_ord_y[self.sites_tot_x[2*j+2]] = \
coor_ord_y[self.sites_tot_x[2*j]] + hop_plt_y[self.bonds_tot_y[j]]
# if the number of AB sites along x increases
l = 2
for j in range(self.sites_cut, sites_dif+2):
for i in range(1+self.sites_tot_x[2*j], self.sites_tot_x[2*j+2]):
coor_ord_x[i] = coor_ord_x[i-1] + self.dx*hop_ord_x[k]
coor_ord_y[i] = coor_ord_y[i-1] + (-1)**l*self.dy*hop_ord_x[k]
k += 1
l += 1
if j < sites_dif+1:
coor_ord_x[self.sites_tot_x[2*j+2]] = coor_ord_x[1+self.sites_tot_x[2*j]]
coor_ord_y[self.sites_tot_x[2*j+2]] = coor_ord_y[1+self.sites_tot_x[2*j]] \
+ hop_plt_y[self.bonds_tot_y[j]]
l = 2
ind_a = -np.ones(self.sites, int)
ind_b = -np.ones(self.sites, int)
for j in range(0, self.sites_cut):
ind_a[self.sites_tot_x[2*j]: self.sites_tot_x[2*j+1]] = \
arange(self.sites_tot_x[2*j]+1, self.sites_tot_x[2*j+2], 2)
for j in range(0, self.sites_cut):
ind_b[self.sites_tot_x[2*j+1]: self.sites_tot_x[2*j+2]] = \
arange(self.sites_tot_x[2*j], self.sites_tot_x[2*j+2], 2)
for j in range(self.sites_cut, sites_dif+2):
ind_a[self.sites_tot_x[2*j]: self.sites_tot_x[2*j+1]] = \
arange(self.sites_tot_x[2*j], self.sites_tot_x[2*j+2], 2)
for j in range(self.sites_cut, sites_dif+2):
ind_b[self.sites_tot_x[2*j+1]: self.sites_tot_x[2*j+2]] = \
arange(self.sites_tot_x[2*j]+1, self.sites_tot_x[2*j+2], 2)
find_ind_b = np.flatnonzero(ind_a==-1)
find_ind_a = np.flatnonzero(ind_b==-1)
ind_b = np.delete(ind_b, find_ind_a)
ind_a = np.delete(ind_a, find_ind_b)
self.coor_hop_x = np.zeros(self.sites)
self.coor_hop_x[self.ind_a] = coor_ord_x[ind_a]
self.coor_hop_x[self.ind_b] = coor_ord_x[ind_b]
self.coor_hop_y = np.zeros(self.sites)
self.coor_hop_y[self.ind_a] = coor_ord_y[ind_a]
self.coor_hop_y[self.ind_b] = coor_ord_y[ind_b]
center_x = sum(self.coor_hop_x)/self.sites
center_y = sum(self.coor_hop_y)/self.sites
self.coor_hop_x -= center_x
self.coor_hop_y -= center_y
[docs]class tri_dislocation(graphene_flake):
def __init__(self, sites_base, sites_cut=0, t_1=1., t_2=1., t_3=1., ext='png'):
graphene_flake.__init__(self, t_1, t_2, t_3, ext=ext)
self.sites_base = sites_base
self.sites_cut = sites_cut
len_x = 2*sites_base+3 # (+4+3)
len_init_x = 2*sites_base
self.sites_x = np.zeros(len_x, int)
self.bonds_x = np.zeros(len_x, int)
self.bonds_tot_x = np.zeros(len_x, int)
self.sites_tot_x = np.zeros(len_x, int)
self.sites_x[1: 2*sites_base: 2] = range(1, sites_base+1)
self.sites_x[2: 2*sites_base: 2] = range(1, sites_base)
self.sites_x[2*sites_base] = sites_base
self.sites_x[2*sites_base+1] = sites_base
self.sites_x[2*sites_base+2] = sites_base-1
for i in range(1, len_x):
self.sites_tot_x[i] = np.sum(self.sites_x[1:i+1])
self.bonds_y = np.arange(sites_base)
self.bonds_tot_y = np.array([np.sum(self.bonds_y[0:i+1]) for i in range(sites_base)])
print('Y BONDS')
print(self.bonds_y)
print(self.bonds_tot_y)
print('X BONDS')
self.bonds_x[1: 2*sites_base-2: 2] = np.arange(1, sites_base)
self.bonds_x[2: 2*sites_base-1: 2] = np.arange(1, sites_base)
self.bonds_x[2*sites_base-1] = sites_base
self.bonds_x[2*sites_base] = sites_base
self.bonds_x[2*sites_base+1] = self.bonds_x[2*sites_base-3]
self.bonds_x[2*sites_base+2] = self.bonds_x[2*sites_base-2]
print(self.bonds_x)
self.bonds_tot_x = np.array([np.sum(self.bonds_x[0:i+1]) for i in range(2*sites_base-1)])
print(self.bonds_tot_x)
self.sites = self.sites_tot_x[-1]
self.sites_init = self.sites_tot_x[2*sites_base-1]
self.coor_x = np.zeros(self.sites)
self.coor_y = np.zeros(self.sites)
y = 0
for i in range(len_init_x-1):
self.coor_x[self.sites_tot_x[i]: self.sites_tot_x[i+1]] = \
self.dx*arange(1-self.sites_x[i+1], self.sites_x[i+1]+1./2, 2.)
self.coor_y[self.sites_tot_x[i]: self.sites_tot_x[i+1]] = y
y += .5
if i%2 == 0:
y += .5
self.sites_disloc = self.sites_init+sites_base
# dislocation
self.coor_y[self.sites_init: self.sites_disloc] = y # - 1 + 1
self.coor_x[self.sites_init: self.sites_disloc: 2] = \
+ self.dx*(1-self.sites_x[2*sites_base-1])+2*self.dx*np.arange(0, sites_base, 2)+0.5
self.coor_x[1+self.sites_init: self.sites_disloc: 2] = \
self.dx*(1-self.sites_x[2*sites_base-1])+2*self.dx*np.arange(0, sites_base, 2)+2*self.dx-0.5
y += self.dx
for i in range(len_init_x, len_x-1):
self.coor_x[self.sites_tot_x[i]: self.sites_tot_x[i+1]] = \
self.dx*arange(1-self.sites_x[i+1], self.sites_x[i+1]+1./2, 2.)
self.coor_y[self.sites_tot_x[i]: self.sites_tot_x[i+1]] = y
y += .5
if i%2 == 0:
y += .5
self.rows_a = len(self.sites_tot_x[0::2])-1
self.rows_b = len(self.sites_tot_x[1::2])
#self.coor_y *= -1.
center_x = sum(self.coor_x)/len(self.coor_x)
self.coor_x -= center_x
center_y = sum(self.coor_y)/len(self.coor_y)
self.coor_y -= center_y
self.get_indices()
fig, ax = plt.subplots()
fig.canvas.set_window_title('LATTICE')
ax.set_xlim([np.min(self.coor_x)-2, np.max(self.coor_x)+2])
ax.set_ylim([np.min(self.coor_y)-2, np.max(self.coor_y)+2])
ax.set_aspect('equal')
plt.plot(self.coor_x[self.ind_a], self.coor_y[self.ind_a], 'bo', markersize=8)
plt.plot(self.coor_x[self.ind_b], self.coor_y[self.ind_b], 'ro', markersize=8)
[docs] def get_eig_lin_strain(self, strain=0, eig_vec=False):
self.strain = strain
# hoppings t_1 (along the y direction)
ind_ham_y_a = np.zeros(self.bonds_tot_y[-1])
ind_ham_y_b = np.zeros(self.bonds_tot_y[-1])
self.hop_y = np.zeros(self.bonds_tot_y[-1])
for i in range(0, self.sites_base-1):
ym_1 = 0.5*(self.coor_y[self.sites_tot_x[2*i]]
+ self.coor_y[self.sites_tot_x[2*i+1]])
print(i, 'i')
ind_ham_y_a[self.bonds_tot_y[i]: self.bonds_tot_y[i+1]] = \
arange(self.sites_tot_x[2*i], self.sites_tot_x[2*i+1])
ind_ham_y_b[self.bonds_tot_y[i]: self.bonds_tot_y[i+1]] = \
arange(self.sites_tot_x[2*i+1], self.sites_tot_x[2*i+2])
self.hop_y[self.bonds_tot_y[i]:self.bonds_tot_y[i+1]] = \
self.t_1*(1.-0.5*strain*ym_1)
# hoppings t_2 and t_3 (along the x direction)
len_x = len(self.sites_tot_x)//2-1
ind_ham_x_a = np.zeros(self.bonds_tot_x[-1])
ind_ham_x_b = np.zeros(self.bonds_tot_x[-1])
self.hop_x = np.zeros(self.bonds_tot_x[-1])
for i in range(1, self.sites_base):
xm_2 = 0.5*(self.coor_x[1+self.sites_tot_x[2*i]: self.sites_tot_x[2*i+1]]
+ self.coor_x[self.sites_tot_x[2*i-1]: self.sites_tot_x[2*i]])
xm_3 = 0.5*(self.coor_x[self.sites_tot_x[2*i]: self.sites_tot_x[2*i+1]-1]
+ self.coor_x[self.sites_tot_x[2*i-1]: self.sites_tot_x[2*i]])
ym = 0.5*(self.coor_y[self.sites_tot_x[2*i-1]]+self.coor_y[self.sites_tot_x[2*i]])
sc_t2 = -self.dx*xm_2-0.5*ym
sc_t3 = self.dx*xm_3-0.5*ym
ind_ham_x_a[self.bonds_tot_x[2*i-2]: self.bonds_tot_x[2*i-1]] = \
arange(self.sites_tot_x[2*i], self.sites_tot_x[2*i+1]-1)
ind_ham_x_b[self.bonds_tot_x[2*i-2]: self.bonds_tot_x[2*i-1]] = \
arange(self.sites_tot_x[2*i-1], self.sites_tot_x[2*i])
self.hop_x[self.bonds_tot_x[2*i-2]: self.bonds_tot_x[2*i-1]] = \
self.t_3*(1.-0.5*strain*sc_t3)
ind_ham_x_a[self.bonds_tot_x[2*i-1]: self.bonds_tot_x[2*i]] = \
arange(1+self.sites_tot_x[2*i], self.sites_tot_x[2*i+1])
ind_ham_x_b[self.bonds_tot_x[2*i-1]: self.bonds_tot_x[2*i]] = \
arange(self.sites_tot_x[2*i-1], self.sites_tot_x[2*i])
self.hop_x[self.bonds_tot_x[2*i-1]: self.bonds_tot_x[2*i]] = \
self.t_2*(1.-0.5*strain*sc_t2)
iA = np.concatenate([ind_ham_x_a, ind_ham_x_b, ind_ham_y_a, ind_ham_y_b])
iB = np.concatenate([ind_ham_x_b, ind_ham_x_a, ind_ham_y_b, ind_ham_y_a])
hop = np.concatenate([self.hop_x, self.hop_x, self.hop_y, self.hop_y])
ham = sparse.csr_matrix((hop, (iA, iB)), shape=(self.sites_init, self.sites_init))
print(ham.todense())
print('solve the eigenvalue problem')
if eig_vec:
self.eig_ener, self.eig_vec = la.linalg.eigh(ham.toarray())
#self.eig_pola = \
#np.array([np.dot(self.eig_vec[self.ind_a, i], self.eig_vec[self.ind_a, i]) for i in range(self.sites)])
else:
self.eig_ener = la.linalg.eigvalsh(ham.toarray())
[docs] def get_lin_strain_lims(self):
"""
Get the maximal negative and positive strain keeping the hoppings positive.
"""
y_m = 0.5*(self.coor_y[0]+self.coor_y[1])
strain_max_pos = +2./y_m - 1e-4
print(strain_max_pos)
i = self.sites_base
y_m = 0.5*(self.coor_y[self.sites_tot_x[2*i-1]]+self.coor_y[self.sites_tot_x[2*i]])
print(self.coor_y[self.sites_tot_x[2*i-1]], self.coor_y[self.sites_tot_x[2*i]])
strain_max_neg = +2./y_m + 1e-4
print('\nMaximal strain values keeping positive hoppings are between:')
print(' ', strain_max_neg, ' ', strain_max_pos, '\n')
return strain_max_neg, strain_max_pos
[docs]def hop_to_real_mw(hop):
"""
Get the sites position in real space from the hoppings.
"""
a = 4.5051e9 # GHz
ga = 0.8287e3 # mm^{−1}
de = 0.00584e9 # GHz
return -2./ga*np.log((hop-de)/a)
# Structured arrays (aka Record array)
# np.recarray((2,), dtype=[('x', int), ('y', float), ('z', int)])
'''
# A site vacancy (not at the center of mass)
flake = tri_to_hexa_zz(14)
strain_max_neg, strain_max_pos = flake.get_lin_strain_lims()
flake.get_eig_lin_strain(strain=strain_max_pos, eig_vec=True,
vacancy=True, row_vacancy=10)
flake.plt_lattice(plt_bonds=False, save=True)
flake.plt_spec(ener_lim=[-2.2, 2.2], nbr_bins=71, save=True)
flake.get_states_selec(e_min=-0.05, e_max=0.05)
flake.plt_states_selec(coef=100000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=True)
flake.get_states_selec(e_min=0.6, e_max=0.9)
flake.plt_states_selec(coef=10000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=True)
flake.get_states_selec(e_min=0.9, e_max=1.3)
flake.plt_states_selec(coef=100000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=True)
flake.get_states_selec(e_min=1.3, e_max=1.5)
flake.plt_states_selec(coef=10000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=True)
flake.get_states_selec(e_min=1.5, e_max=1.7)
flake.plt_states_selec(coef=100000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=True)
flake.get_states_selec(e_min=1.7, e_max=1.9)
flake.plt_states_selec(coef=10000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=True)
flake.get_states_selec(e_min=1.9, e_max=2.)
flake.plt_states_selec(coef=100000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=True)
flake.get_states_selec(e_min=2., e_max=2.1)
flake.plt_states_selec(coef=10000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=True)
'''
# B site vacancy centered #i = 9
'''
flake = tri_to_hexa_zz(12)
strain_max_neg, strain_max_pos = flake.get_lin_strain_lims()
flake.get_eig_lin_strain(strain=strain_max_pos, eig_vec=True,
vacancy=True, row_vacancy=9)
flake.plt_lattice(plt_bonds=False, save=True)
flake.plt_spec(ener_lim=[-2.5, 2.5], nbr_bins=51, save=True)
flake.get_states_selec(e_min=-0.05, e_max=0.05)
flake.plt_states_selec(coef=100000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=True)
flake.get_states_selec(e_min=0.9, e_max=1.4)
flake.plt_states_selec(coef=100000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=True)
flake.get_states_selec(e_min=1.4, e_max=1.6)
flake.plt_states_selec(coef=10000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=True)
flake.get_states_selec(e_min=1.6, e_max=1.8)
flake.plt_states_selec(coef=100000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=True)
flake.get_states_selec(e_min=1.8, e_max=2.05)
flake.plt_states_selec(coef=10000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=True)
flake.get_states_selec(e_min=2.05, e_max=2.2)
flake.plt_states_selec(coef=100000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=True)
flake.get_states_selec(e_min=2.2, e_max=2.3)
flake.plt_states_selec(coef=10000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=True)
'''
'''
# A site vacancy (not at the center of mass)
flake = tri_to_hexa_zz(13) #i = 8
strain_max_neg, strain_max_pos = flake.get_lin_strain_lims()
flake.get_eig_lin_strain(strain=strain_max_pos, eig_vec=True, vacancy=True, row_vacancy=8)
flake.plt_lattice(plt_bonds=False, save=False)
flake.plt_spec(ener_lim=[-2., 2.], nbr_bins=51)
flake.get_states_selec(e_min=-0.05, e_max=0.05)
flake.plt_states_selec(coef=100000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=False)
flake.get_states_selec(e_min=0.8, e_max=0.9)
flake.plt_states_selec(coef=10000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=False)
flake.get_states_selec(e_min=0.9, e_max=1.3)
flake.plt_states_selec(coef=100000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=False)
flake.get_states_selec(e_min=1.3, e_max=1.6)
flake.plt_states_selec(coef=10000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=False)
flake.get_states_selec(e_min=1.6, e_max=1.7)
flake.plt_states_selec(coef=100000., fontsize=20, plt_sites=True,
plt_bonds=False, plt_axis=False, save=False)
'''
plt.show()