Kwant-discuss

kwant-discuss@python.org

918 discussions

Dispersion Plot of SSH Model with Linearly Varying Chemical Potential
by Subhadeep Chakraborty 27 Jun '22

27 Jun '22

Hi Kwant Users ! I was trying to plot the band structure of the ssh model with linearly varying chemical potential (\mu*x). But it is showing error that "*Expecting an arraylike object or a scalar*" while adding this potential term. Can anyone please tell where I am making mistake and how to rectify it ? _____________________________________ Code is attached below :- import kwant import numpy as np import math ## import matplotlib #matplotlib.use('Agg') ## # For plotting from matplotlib import pyplot as plt # For matrix support import tinyarray I2 = tinyarray.array([[1,0], [0,1]]) Sx = tinyarray.array([[0,1], [1,0]]) Sy = tinyarray.array([[0,-1j], [1j,0]]) Sz = tinyarray.array([[1,0], [0,-1]]) Sp = tinyarray.array([[0,2], [0,0]]) Sm = tinyarray.array([[0,0], [2,0]]) def create_system(length=10): t1 = 0.5 ; t2 = 1.0 ; def Onsite(site): (x, ) = site.pos potential = x return x # system building lat = kwant.lattice.square(a=1, norbs=2) syst = kwant.Builder() # central scattering region syst[(lat(x, 0) for x in range(length))] = t1*Sx syst[lat.neighbors()] = (t2*Sx - 1j*t2*Sy)/2 # add leads sym = kwant.TranslationalSymmetry((-1, 0)) lead_left = kwant.Builder(sym) lead_left[lat(0, 0)] = t1*Sx + Onsite*I2 lead_left[lat.neighbors()] = (t2*Sx - 1j*t2*Sy)/2 syst.attach_lead(lead_left) syst.attach_lead(lead_left.reversed()) return syst def main(): # parameters # create system syst = create_system(length=10).finalized() # plot the system and dispersion kwant.plot(syst) kwant.plotter.bands(syst.leads[0], show=False) plt.plot([-np.pi, np.pi], [chemical_potential, chemical_potential], 'k--') plt.show() lead = syst.leads[0] bands = kwant.physics.Bands(lead) momenta = np.linspace(-2*np.pi, 2*np.pi, 401) energies = [bands(k) for k in momenta] plt.plot(momenta, energies,linestyle='--') plt.xlim(-2*np.pi, 2*np.pi) #pyplot.ylim(-5, 5) plt.xlabel("$k$") plt.ylabel("Energy") #pyplot.axhline(y=0.0, color='k', linestyle='-') plt.tight_layout() plt.grid() plt.show() plt.close() if __name__ == '__main__': main() ________________________________________________ Thanks in advance, *_________________________________* *Subhadeep Chakraborty* *_________________________________*

2 1

Plot Band-Structure of SSH Model with Linearly Varying Chemical Potential Using Kwant.continuum
by Subhadeep Chakraborty 27 Jun '22

27 Jun '22

Hi Kwant Users ! I was trying to plot the band structure of the SSH model using kwant.continuum. I could do it when the hopping terms are constant. Now let say the hoppings are linearly varying along x direction (t*x) or say I want to add a chemical potential like (\mu*x). In that case, how should I introduce this potential in the code ? I have attached my code below. It will be a great help if anyone kindly tells me where I need to modify the code and how to insert this (\mu*x term) in the code. The code is attached herewith. _______________________________________ import kwant import kwant.continuum import scipy.sparse.linalg import scipy.linalg import numpy as np # For plotting import matplotlib as mpl from matplotlib import pyplot as plt def ssh(): hamiltonian = "m*sigma_x + alpha * k_x * sigma_y" def plot(ax, a=1): params = dict(m=.5, alpha=.5) h_k = kwant.continuum.lambdify(hamiltonian, locals=params) k_cont = np.linspace(-1, 1, 201) e_cont = [scipy.linalg.eigvalsh(h_k(k_x=ki)) for ki in k_cont] template = kwant.continuum.discretize(hamiltonian, grid=a) syst = kwant.wraparound.wraparound(template).finalized() def h_k(k_x): p = dict(k_x=k_x, **params) return syst.hamiltonian_submatrix(params=p) ax.plot(k_cont, e_cont, 'r-') ax.plot([], [], 'r-', label='continuum') ax.set_xlim(-1, 1) ax.set_ylim(-2, 2) ax.set_title('a={}'.format(a)) ax.set_xlabel('Momentum k_x') ax.set_ylabel('energy E') ax.grid() ax.legend() _, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12, 4)) plot(ax1, a=1) plot(ax2, a=.25) plt.show() def main(): ssh() if __name__ == '__main__': main() _______________________________________________________ *_Thanks in advance,*

3 2

Removing a certain hopping in the Kwant.plot
by Sayandip Dhara 03 Jun '22

03 Jun '22

Hello, As I understand that when I use "kwant.plot" to generate a lattice it automatically connects with a line whenever it sees a hopping between two sites. Now in my case, I wanted to connect two edges of a system by hopping to implement periodic boundary conditions. However, now when I plot the lattice I can see the long black lines from one edge to another since they are connected by a hopping. So, when I do a ldos I can still see those lines in my plot. So is there a way to remove the lines corresponding to the hopping between two edges in the plot?

3 3

Removing a certain hopping in the Kwant.plot
by Sayandip Dhara 03 Jun '22

03 Jun '22

1 0

How to add random but opposite onsite potential on graphene?
by wilson2048＠outlook.com 26 May '22

26 May '22

Dear community, I want to add random potentials on graphene. Within each unit cell, I want A site have the oppsite value to B site. That is, random across unit cells, but symmetric respect to zero within every unit cell. How to achieve this? I have tried to do so, but cannot find a way. Below is my code: ``` import numpy as np import kwant low = 0 high = 1 def make_syst(a=1): syst = kwant.Builder() lat = kwant.lattice.honeycomb(a, norbs=1, name=['a', 'b']) r = 10 def circle(pos): x, y = pos return x ** 2 + y ** 2 < r ** 2 def onsite(site, low, high): return np.random.uniform(low, high) # first attempt using function, not working, also cannot gurantee a, b are truly opposite in value # syst[lat.a.shape(circle, (0, 0))] = onsite # syst[lat.b.shape(circle, (0, 0))] = -1*onsite # wrong # second attempt, still not right m = np.random.uniform(low, high) syst[lat.a.shape(circle, (0, 0))] = m syst[lat.b.shape(circle, (0, 0))] = -m syst[lat.neighbors()] = 1 return syst.finalized() fsyst = make_syst() # kwant.plot(fsyst) H = fsyst.hamiltonian_submatrix(params=dict(high=high, low=low)) Es = np.linalg.eigvalsh(H) plt.plot(Es, '.') plt.show() ``` Thanks for help!

3 3

Kwant-Plot discretize potential with different colours
by Lorenzo BAGNASACCO 17 May '22

17 May '22

Dear Kwant community, I am a physics student and I am approaching to Kwant to solve simple problems of scattering. For example, I created a 2D discrete potential. The potential has two different values: "V0" and "0". When I plot the system, using the command "kwant.plot(syst)", I obtain the plot with only light-blue colour. There is a command to plot the discrete potential with two (or more) different colours (in my case one colour for "V0", and another one for "0")? Yours Sincerely Lorenzo Bagnasacco

2 1

Problems in getting the current
by mondal18＠iitg.ac.in 12 May '22

12 May '22

Dear developers, I was trying to get the plot of currents in a graphene ribbon. I am getting the error "Number of orbitals not defined", which is done by specifying the norbs. However, in my code, I am unable to specify the norbs. Another thing is that the interface between the right lead (and left leads) and the system is not parallel to the y axis. What are the possible ways to solve such problems. The code I am using is the following. #################################################### from math import sqrt import random import itertools as it import tinyarray as ta import numpy as np import matplotlib.pyplot as plt import kwant class Honeycomb(kwant.lattice.Polyatomic): """Honeycomb lattice with methods for dealing with hexagons""" def __init__(self, name=''): prim_vecs = [[0.5, sqrt(3)/2], [1, 0]] # bravais lattice vectors # offset the lattice so that it is symmetric around x and y axes basis_vecs = [[-0.5, -1/sqrt(12)], [-0.5, 1/sqrt(12)]] super(Honeycomb, self).__init__(prim_vecs, basis_vecs, name) self.a, self.b = self.sublattices def hexagon(self, tag): """ Get sites belonging to hexagon with the given tag. Returns sites in counter-clockwise order starting from the lower-left site. """ tag = ta.array(tag) # a-sites b-sites deltas = [(0, 0), (0, 0), (1, 0), (0, 1), (0, 1), (-1, 1)] lats = it.cycle(self.sublattices) return (lat(*(tag + delta)) for lat, delta in zip(lats, deltas)) def hexagon_neighbors(self, tag, n=1): """ Get n'th nearest neighbor hoppings within the hexagon with the given tag. """ hex_sites = list(self.hexagon(tag)) return ((hex_sites[(i+n)%6], hex_sites[i%6]) for i in range(6)) def random_placement(builder, lattice, density): """ Randomly selects hexagon tags where adatoms can be placed. This avoids the edge case where adatoms would otherwise be placed on incomplete hexagons at the system boundaries. """ assert 0 <= density <= 1 system_sites = builder.sites() def hexagon_in_system(tag): return all(site in system_sites for site in lattice.hexagon(tag)) # get set of tags for sites in system (i.e. tags from only one sublattice) system_tags = (s.tag for s in system_sites if s.family == lattice.a) # only allow tags which have complete hexagons in the system valid_tags = list(filter(hexagon_in_system, system_tags)) N = int(density * len(valid_tags)) total_hexagons=len(valid_tags) valuef=random.sample(valid_tags, N) return valuef def ribbon(W, L): def shape(pos): return (-L <= pos[0] <= L and -W <= pos[1] <= W) return shape ## Pauli matrices ## sigma_0 = ta.array([[1, 0], [0, 1]]) # identity sigma_x = ta.array([[0, 1], [1, 0]]) sigma_y = ta.array([[0, -1j], [1j, 0]]) sigma_z = ta.array([[1, 0], [0, -1]]) ## Hamiltonian value functions ## def onsite_potential(site, params): return params['ep'] * sigma_0 def potential_shift(site, params): return params['mu'] * sigma_0 def kinetic(site_i, site_j, params): return -params['gamma'] * sigma_0 def rashba(site_i, site_j, params): d_ij = site_j.pos - site_i.pos rashba = 1j * params['V_R'] * (sigma_x * d_ij[1] - sigma_y * d_ij[0]) return rashba + kinetic(site_i, site_j, params) def spin_orbit(site_i, site_j, params): so = 1j * params['V_I'] * sigma_z return so lat = Honeycomb() A, B = lat.sublattices pv1, pv2 = lat.prim_vecs ysym = kwant.TranslationalSymmetry(pv2 - 2*pv1) # lattice symmetry in -y direction xsym = kwant.TranslationalSymmetry(-pv2) # lattice symmetry in -x direction # adatom lattice, for visualization only adatom_lat = kwant.lattice.Monatomic(lat.prim_vecs, name='adatom') def site_size(s): return 0.25 if s.family in lat.sublattices else 0.3 def site_color(s): if s.family==lat.sublattices[0]: return 'w' elif s.family==lat.sublattices[1]: return 'k' else: return 'g' def create_lead_h(W, sym1, axis=(0, 0)): lead = kwant.Builder(sym1) lead[lat.wire(axis, W)] = 0. * sigma_0 lead[lat.neighbors(1)] = kinetic return lead def create_ribbon(W, L, adatom_density=0.2, show_adatoms=False): ## scattering region ## sys = kwant.Builder() sys[lat.shape(ribbon(W, L), (0, 0))] = onsite_potential sys[lat.neighbors(1)] = kinetic adatoms = random_placement(sys, lat, adatom_density) sys[(lat.hexagon(a) for a in adatoms)] = potential_shift sys[(lat.hexagon_neighbors(a, 1) for a in adatoms)] = rashba sys[(lat.hexagon_neighbors(a, 2) for a in adatoms)] = spin_orbit if show_adatoms: # no hoppings are added so these won't affect the dynamics; purely cosmetic sys[(adatom_lat(*a) for a in adatoms)] = 0 ## leads ## leads = [create_lead_h(W, xsym)] leads += [lead.reversed() for lead in leads] # right lead for lead in leads: sys.attach_lead(lead) return sys def plot_currents(syst, currents): fig, axes = plt.subplots(1, len(currents)) if not hasattr(axes, '__len__'): axes = (axes,) for ax, (title, current) in zip(axes, currents): kwant.plotter.current(syst, current, ax=ax, colorbar=False) ax.set_title(title) plt.show() def currents(sys, params): syst = sys.finalized() wf = kwant.wave_function(syst, energy=-1, args=(params,)) psi = wf(0)[0] # even (odd) indices correspond to spin up (down) up, down = psi[::2], psi[1::2] print(psi.shape) density = np.abs(up)**2 + np.abs(down)**2 # spin down components have a minus sign spin_z = np.abs(up)**2 - np.abs(down)**2 # spin down components have a minus sign spin_y = 1j * (down.conjugate() * up - up.conjugate() * down) rho = kwant.operator.Density(syst) rho_sz = kwant.operator.Density(syst, sigma_z) rho_sy = kwant.operator.Density(syst, sigma_y) # calculate the expectation values of the operators with 'psi' density = rho(psi) spin_z = rho_sz(psi) spin_y = rho_sy(psi) plot_densities(syst, [ ('$σ_0$', density), ('$σ_z$', spin_z), ('$σ_y$', spin_y), ]) J_0 = kwant.operator.Current(syst) J_z = kwant.operator.Current(syst, sigma_z) J_y = kwant.operator.Current(syst, sigma_y) # calculate the expectation values of the operators with 'psi' current = J_0(psi) spin_z_current = J_z(psi) spin_y_current = J_y(psi) plot_currents(syst, [ ('$J_{σ_0}$', current), ('$J_{σ_z}$', spin_z_current), ('$J_{σ_y}$', spin_y_current), ]) if __name__ == '__main__': params = dict(gamma=1., ep=0, mu=0., V_I=0.007, V_R=0.0165) # In adatoms W=11 L=32 adatom_density=0.1 sys = create_ribbon(W, L, adatom_density, show_adatoms=True) # fig = plt.figure() ax = plt.subplot() ax = kwant.plot(sys, site_color=site_color, site_lw=0.05, site_size=site_size,lead_site_lw=0, dpi = 150) # ax.savefig('system2.png', bbox_inches = 'tight', dpi = 200) currents(sys, params) ################################################### I shall look forward to your kind reply. Thanking you. Sincerely yours, Sayan Mondal, Research scholar, IIT Guwahati, India

2 1

Help: Inputting Random Potential on Graphene
by Blackwell, Raymond 12 May '22

12 May '22

Hello, My name is Raymond Blackwell and I am attempting to use Kwant to simulate some experimental data. Our data involves dopant atoms intercalated under graphene, so I would like to simulate the LDOS but with a potential that randomly varies at several points within the system. I am trying to modify the code for the kernel_polynominal_method, but I cannot seem to get a random potential. I have used a similar method to randomly select a few points and modify the hopping parameter, but I am struggling with the seemingly more simple case where the potential varies. The relevant portion of my code is below. def make_syst(r=30, t=-1, a=1): syst = kwant.Builder() lat = kwant.lattice.honeycomb(a, norbs=1) def circle(pos): x, y = pos return x ** 2 + y ** 2 < r ** 2 sites2 = list(syst.sites()) Random_sites= random.choices(sites2, k=10) def onsite(site, salt): return 0.5 * kwant.digest.gauss(repr(Random_sites), salt) syst[lat.shape(circle, (0, 0))] = onsite(bytes(3),bytes(3)) Any help is greatly appreciated! Best, -Raymond

2 3

Re: Problem in calculation of Current density
by Ajit Kumar Sahu 11 Apr '22

11 Apr '22

1 0

current-voltage diagram
by loojoo026＠gmail.com 08 Apr '22

08 Apr '22

dear all I have done my project using Kwant and I have drawn a conductance-energy diagram for the system. The question I have is how can I get the current-voltage diagram? thanks Leo

4 4

Jump to page:

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

Kwant-discuss