Mailman 3 October 2021 - Kwant-discuss

Interpolation is slow
by Mani Chandra 30 Oct '21

30 Oct '21

Hello, We're dealing with a 5000 x 5000 site system, and the time to interpolate the densities and currents using kwant.plotter.interpolate_density and kwant.plotter.interpolate_current are by far the rate limiting steps in the calculation: System assembly : 15 mins Computing scattering states: 3 hrs Interpolating density and current (for use in plotting): 6 hrs! Is there any way to speed up the interpolation routines? I see that the present implementation in plotter.py involves a for loop (in python) over each site/hopping. Perhaps a vectorized implementation may help? Thanks, Mani

2 8

Conductivity calculation using the kpm method for a spinful system
by mondal18＠iitg.ac.in 29 Oct '21

29 Oct '21

Dear developer, I was trying to calculate the conductivity of honeycomb lattice, where I have included the Haldane term along with the spin-orbit Rashba term (Kane-Mele model). I was following the tutorial provided in chapter 2.9 (https://kwant-project.org/doc/1/tutorial/kpm). The code shows an error when I include the spin in the system. My code is attached below, from matplotlib import pyplot as plt import tinyarray as tn import kwant import numpy as np t = 1 t2 = 0.3 sgm_0=tn.array([[1,0],[0,1]]) sgm_x=tn.array([[0,1],[1,0]]) sgm_y=tn.array([[0,-1j],[1j,0]]) sgm_z=tn.array([[1,0],[0,-1]]) def make_syst_topo(l=20, a=1, m=0): alpha = 0.01 syst = kwant.Builder() lat = kwant.lattice.honeycomb(a, norbs=1, name=['a', 'b']) a = lat.a b = lat.b def rectangle(pos): x, y = pos return - l<= x <= +l and - l<= y <= +l def rashba(site_i,site_j): d_ij = site_i.pos - site_j.pos return -1j*alpha*(sgm_x*d_ij[1]-sgm_y*d_ij[0]) def nn_hop(site_i, site_j): return rashba(site_i, site_j) - t * sgm_0 syst[lat.shape(rectangle, (0, 0))] = 0.0*sgm_0 syst[kwant.builder.HoppingKind((0, 0), a, b)] = nn_hop syst[kwant.builder.HoppingKind((0, 1), a, b)] = nn_hop syst[kwant.builder.HoppingKind((-1, 1), a, b)] = nn_hop syst[lat.a.neighbors()] = 1j * t2 * sgm_z syst[lat.b.neighbors()] = -1j * t2*sgm_z syst.eradicate_dangling() return lat, syst.finalized() def plot_system(syst): kwant.plot(syst,site_size=0.1, hop_lw=0.01) # Plot fill density of states plus curves on the same axes. def plot_dos_and_curves(dos, labels_to_data): plt.figure(figsize=(10,8)) plt.fill_between(dos[0], dos[1], label="DoS [a.u.]", alpha=0.5, color='gray') for label, (x, y) in labels_to_data: plt.plot(x, y, label=label, linewidth=2) plt.legend(loc='upper center', framealpha=0.5, fontsize=12) plt.xlabel(r'$E/t$'); # plt.ylabel('ylabel') plt.show() plt.clf() def conductivity_example(): # construct the Haldane model lat, fsyst = make_syst_topo() # find 'A' and 'B' sites in the unit cell at the center of the disk where = lambda s: np.linalg.norm(s.pos) < 1 # component 'xx' s_factory = kwant.kpm.LocalVectors(fsyst, where) cond_xx = kwant.kpm.conductivity(fsyst, alpha='x', beta='x', mean=True, num_vectors=None, vector_factory=s_factory) # component 'xx' s_factory = kwant.kpm.LocalVectors(fsyst, where) cond_yy = kwant.kpm.conductivity(fsyst, alpha='y', beta='y', mean=True, num_vectors=None, vector_factory=s_factory) # component 'xy' s_factory = kwant.kpm.LocalVectors(fsyst, where) cond_xy = kwant.kpm.conductivity(fsyst, alpha='x', beta='y', mean=True, num_vectors=None, vector_factory=s_factory) energies_x = cond_xx.energies energies_y = cond_yy.energies cond_array_xx = np.array([cond_xx(e, temperature=0.00) for e in energies_x]) cond_array_xy = np.array([cond_xy(e, temperature=0.00) for e in energies_x]) cond_array_yy = np.array([cond_yy(e, temperature=0.00) for e in energies_y]) # area of the unit cell per site area_per_site = np.abs(np.cross(*lat.prim_vecs)) / len(lat.sublattices) cond_array_xx /= area_per_site cond_array_xy /= area_per_site cond_array_yy /= area_per_site # ldos s_factory = kwant.kpm.LocalVectors(fsyst, where) spectrum = kwant.kpm.SpectralDensity(fsyst, num_vectors=None, vector_factory=s_factory) plot_dos_and_curves( (spectrum.energies, spectrum.densities * 8), [ (r'$\sigma_{xx} / 4$', (energies_x, cond_array_xx / 4)), (r'$\sigma_{xy}$', (energies_x, cond_array_xy)), (r'$\sigma_{yy} / 4$', (energies_y, cond_array_yy / 4))] ) def main(): fsyst = make_syst_topo()[1] plot_system(fsyst) conductivity_example() if __name__ == '__main__': main() I shall look forward to your kind reply.

1 0

How to add leads oblique to a square lattice?
by Jerry xhm 13 Oct '21

13 Oct '21

Hello Kwant community, Given a square lattice Hamiltonian, I am trying to put a rectangular sample along an oblique direction to the x,y axes of the underlying square lattice, i.e., the rectangular sides are along new rotated X,Y axes. One can fix the site at the origin (0,0) in the rotation. I have done this successfully using commensurate rotation angles in the following code. The comments in the code explain in detail. The sample looks like this (https://www.dropbox.com/s/2rd9bjbd7ataosy/lattice_fig.pdf?dl=0). Then I would like to attach leads along the new X-axis (just normally connect to the rectangular sample that doesn't need adding extra sites). I put the translation symmetry rotated to X and confine the lead shape in Y. I expected the lead would just connect to the sample with the interface in Y direction. However, many extra sites are added and the interface is rendered in the old y direction, although the lead does extend in X. The code below attaches only one lead in order to contrast the interface and the opposite unaffected sample edge. The system looks like this (https://www.dropbox.com/s/3wj6r2ws5ev20p5/lattice_fig0.pdf?dl=0). So the question is how to attach leads in the expected way. Thanks!

3 8

Graphene Sparse vs Dense - Energy calculation values dependent on graphene length
by Henry Axt 06 Oct '21

06 Oct '21

This is related to my other post, but as I see them as different (maybe not?) problems, I will post a second thread. The situation is, again, a graphene strip with periodic boundary conditions in one direction and open boundaries in the other. Using the sparse solver from linalg, I have solved for energies and gotten reasonable results for the energies, except for the case when the length (along the axis of the open boundary condition) is metallic. Attached is a code that shows it for length = 17, which is a metallic length for graphene (has energies at zero), but it does not produce messy energy calculation for when it is not a metallic length, i.e. something like length = 16. """ Graphene wire with periodic boundary conditions in the vertical direction. """ import kwant from math import pi, sqrt, tanh, cos, ceil, floor, atan, acos, asin from cmath import exp import numpy as np import scipy import scipy.linalg as lina sin_30, cos_30, tan_30 = (1 / 2, sqrt(3) / 2, 1 / sqrt(3)) def create_closed_system(length, width, lattice_spacing, onsite_potential, hopping_parameter, boundary_hopping): padding = 0.5*lattice_spacing*tan_30 def calc_total_length(length): total_length = length N = total_length//lattice_spacing # Number of times a graphene hexagon fits in horizontially fully new_length = N*lattice_spacing + lattice_spacing*0.5 diff = total_length - new_length if diff != 0: length = length - diff total_length = new_length return total_length def calc_width(width,lattice_spacing, padding): stacking_width = lattice_spacing*((tan_30/2)+(1/(2*cos_30))) N = width//stacking_width if N % 2 == 0.0: # Making sure that N is odd. N = N-1 new_width = N*stacking_width + padding width = new_width return width, int(N) def rectangle(pos): x,y = pos if (0 < x <= total_length) and (-padding <= y <= width -padding): return True return False def lead_shape(pos): x, y = pos return 0 - padding <= y <= width def tag_site_calc(x): return int(-1*(x*0.5+0.5)) #Initation of geometrical limits of the lattice of the system total_length = calc_total_length(length) width, N = calc_width(width, lattice_spacing, padding) # The definition of the potential over the entire system def potential_e(site,pot): return pot ### Definig the lattices ### graphene_e = kwant.lattice.honeycomb(a=lattice_spacing,name='e') a, b = graphene_e.sublattices sys = kwant.Builder() # The following functions are required for input in the def onsite_shift_e(site, pot): return potential_e(site,pot) sys[graphene_e.shape(rectangle, (0.5*lattice_spacing, 0))] = onsite_shift_e sys[graphene_e.neighbors()] = -hopping_parameter ### Boundary conditions for scattering region ### for site in sys.sites(): (x,y) = site.tag if float(site.pos[1]) < 0: if str(site.family) == "<Monatomic lattice e1>": sys[b(x,y),a(int(x+tag_site_calc(N)),N)] = -boundary_hopping kwant.plot(sys) return sys def eigen_vectors_and_values(sys,sparse_dense,k): #sparse = 0, dense = 1 if sparse_dense == 0: ham_mat = sys.hamiltonian_submatrix(params=dict(pot=0.0), sparse=True) eigen_val, eigen_vec = scipy.sparse.linalg.eigsh(ham_mat.tocsc(), k=k, sigma=0, return_eigenvectors=True) if sparse_dense == 1: ham_mat = sys.hamiltonian_submatrix(params=dict(pot=0.0), sparse=False) eigen_val, eigen_vec = lina.eigh(ham_mat) #sort the ee and ev in ascending order idx = eigen_val.argsort() eigen_val= eigen_val[idx] eigen_vec = eigen_vec[:,idx] if sparse_dense == 1: eigen_val = eigen_val[int(len(eigen_val)*0.5-k*0.5):] # Remove first part eigen_val = eigen_val[:k] # Take first k-values eigen_vec = eigen_vec[:,int(len(eigen_val)*0.5-k*0.5):] eigen_vec = eigen_vec[:,:k] return eigen_vec, eigen_val def main(): sys = create_closed_system(length = 17.0, width = 20.0, lattice_spacing = 1.0, onsite_potential = 0.0, hopping_parameter = 1.0, boundary_hopping = 1.0) sys = sys.finalized() #sparse = 0, dense = 1 eigen_vec, eigen_val = eigen_vectors_and_values(sys,sparse_dense=0,k=50) eigen_vec2, eigen_val2 = eigen_vectors_and_values(sys,sparse_dense=1,k=50) print("Eigenenergy difference (between sparse and dense): ") for i_ in range(len(eigen_val)): print(eigen_val[i_] - eigen_val2[i_] ) return if __name__ == "__main__": main()

2 2