#!/usr/bin/env python from __future__ import print_function # python3 style print # Chain with alternating site energies and hoppings # Study surface properties of finite chain from pythtb import * import matplotlib as mpl import matplotlib.pyplot as plt # to set up model for given surface energy shift and lambda def set_model(n_cell,en_shift,lmbd): # set parameters of model t=-1.0 # average hopping Delta=-0.4*np.cos(lmbd) # site energy alternation del_t=-0.3*np.sin(lmbd) # bond strength alternation # construct bulk model lat=[[1.0]] orb=[[0.0],[0.5]] bulk_model=tbmodel(1,1,lat,orb) bulk_model.set_onsite([Delta,-Delta]) bulk_model.add_hop(t+del_t, 0, 1, [0]) bulk_model.add_hop(t-del_t, 1, 0, [1]) # cut chain of length n_cell and shift energy on last site finite_model=bulk_model.cut_piece(n_cell,0) finite_model.set_onsite(en_shift,ind_i=2*n_cell-1,mode='add') return finite_model # set Fermi energy and number of cells Ef=0.18 n_cell=20 n_orb=2*n_cell # set number of parameter values to run over n_param=101 # initialize arrays params=np.linspace(0.,1.,n_param) eig_sav=np.zeros((n_orb,n_param),dtype=float) xbar_sav=np.zeros((n_orb,n_param),dtype=float) nocc_sav=np.zeros((n_param),dtype=int) surf_sav=np.zeros((n_param),dtype=float) count=np.zeros((n_orb),dtype=float) # initialize plots mpl.rc('font',size=10) # set global font size fig,ax=plt.subplots(3,2,figsize=(7.,6.), gridspec_kw={'height_ratios':[2,1,1]},sharex="col") # loop over two cases: vary surface site energy, or vary lambda for mycase in ['surface energy','lambda']: if mycase == 'surface energy': (ax0,ax1,ax2)=ax[:,0] # axes for plots in left panels ax0.text(-0.30,0.90,'(a)',size=22.,transform=ax0.transAxes) lmbd=0.15*np.pi*np.ones((n_param),dtype=float) en_shift=-3.0+6.0*params abscissa=en_shift elif mycase == 'lambda': (ax0,ax1,ax2)=ax[:,1] # axes for plots in right panels ax0.text(-0.30,0.90,'(b)',size=22.,transform=ax0.transAxes) lmbd=params*2.*np.pi en_shift=0.2*np.ones((n_param),dtype=float) abscissa=params # loop over parameter values for j in range(n_param): # set up and solve model; store eigenvalues my_model=set_model(n_cell,en_shift[j],lmbd[j]) (eval,evec)=my_model.solve_all(eig_vectors=True) # find occupied states nocc=(eval