Hello, My name is Enrique and first of all, thanks for the developping of such a useful tool. My question is related to that one of the edge potentials. I am trying to reproduce the Haldane model with disorder in a rectangular scattering region. My goal is to introduce edge disorder as well as compute the local current density. I have two problems: a) Regarding the edge disorder, I cannot remove the disorder in the edges with the leads (I saw your recent previous discussion) b) Is it possible already to compute directly from Kwant the local current density? If yes, how? Thanks again, Enrique from __future__ import division # so that 1/2 == 0.5, and not 0 from scipy.sparse import spdiags from math import pi, sqrt, tanh, e import numpy as np import kwant import random import math # For computing eigenvalues import scipy.sparse.linalg as sla # For plotting from matplotlib import pyplot # Define the graphene lattice graphene = kwant.lattice.honeycomb() a, b = graphene.sublattices edge_potential=10 def make_system(L=10, W=10, pot=0.001): t=1 tt= 0.04 * t phi=1*pi/2 edge=1 # 1 if there is edge disorder bulk=0 # 1 if there is bulk disorder nnn_hoppinga=tt*e**(1j*phi) nnn_hoppingb=tt*e**(1j*phi) #### Define the scattering region. #### def rec(pos): x, y = pos return 0 < x < L and 0 < y < W sys = kwant.Builder() # w: width and pot: potential maximum of the p-n junction def potential(site): (x, y) = site.pos d = y * 2 + x * 3 return 0 sys[graphene.shape(rec, (1, 1))] = potential # specify the hoppings of the graphene lattice nn_hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b)) nnn_hoppings_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a)) nnn_hoppings_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b)) sys[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb # Modify the scattering region if edge==1: # change the values of the potential at the edge for site in sys.expand(graphene.shape(rec, (0, 1))): if sys.degree(site)<=9 and abs(site.pos[0])!=L-4: #abs(site.pos[0])!=19 execludes the interfaces with the leads sys[site]=edge_potential elif bulk==1: # change the values of the potential at the edge for site in sys.expand(graphene.shape(rec, (0, 1))): if sys.degree(site)==9 and abs(site.pos[0])!=19: #abs(site.pos[0])!=19 execludes the interfaces with the leads sys[site]=edge_potential #### Define the leads. #### # left lead sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0))) sym0.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)]) sym0.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)]) def lead0_shape(pos): x, y = pos return (0 < y < W) lead0 = kwant.Builder(sym0) lead0[graphene.shape(lead0_shape, (0, 0))] = -pot lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb # right lead sym1 = kwant.TranslationalSymmetry(graphene.vec((1, 0))) sym1.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)]) sym1.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)]) def lead1_shape(pos): x, y = pos return (0 < y < W) lead1 = kwant.Builder(sym1) lead1[graphene.shape(lead1_shape, (0, 0))] = pot lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb sys.attach_lead(lead1) sys.attach_lead(lead0) return sys # Call the main function if the script gets executed (as opposed to imported). if __name__ == '__main__': main()t.show() def main(): sys = make_system() # Check that the system looks as intended. def family_color(site): #print(sys[site]) if sys[site]==edge_potential: return 'green' else: return 'black' def site_size(site): if sys[site]==edge_potential: return 0.35 else:return 0.2 kwant.plot(sys,site_color=family_color,site_size=site_size) # Finalize the system. sys = sys.finalized() # Compute number of neighbors i=sys.sites.index(graphene(5,6)) all_the_neighbors=sys.graph.out_neighbors(i) # Call the main function if the script gets executed (as opposed to imported). # See <http://docs.python.org/library/__main__.html><http://docs.python.org/library/__main__.html>. if __name__ == '__main__': main()
Hi,
a) Regarding the edge disorder, I cannot remove the disorder in the edges with the leads (I saw your recent previous discussion)
I am not sure what you mean by this. From what I can see from your model there is no disorder in your leads: lead0 = kwant.Builder(sym0) lead0[graphene.shape(lead0_shape, (0, 0))] = -pot lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb Also, you do not need to manually enumerate all the hopping kinds like this. Lattices have a method 'neighbors' that returns the hopping kinds for neighboring sites in the lattice. You could reduce the above code to: lead0[graphene.shape(lead0_shape, (0, 0))] = -pot lead0[graphene.neighbors()] = -t lead0[(kind.family_a == a for kind in graphene.neighbors(2))] = -nnn_hoppinga lead0[(kind.family_a == b for kind in graphene.neighbors(2))] = -nnn_hoppingb which would allow you to remove the lines where you explicitly enumerate the hoppings: nn_hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b)) nnn_hoppings_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a)) nnn_hoppings_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b))
Is it possible already to compute directly from Kwant the local current density? If yes, how?
This will be available in Kwant version 1.3. I see that at the start of your script that you have the line 'from __future__ import division' which tells me that you are still using Python 2. Since version 1.2 Kwant supports Python 3 only, and the 1.1.x versions will only receive bug fixes, no new features. For this reason I would recommend using Python 3 unless you have a good reason not to. Coming back to your question, while calculating currents will be available in Kwant 1.3, if your model only has 1 orbital per site, it is probably just as easy to calculate the current manually. The current flowing from site 'j' to site 'i' is written: J_ij = 2 Im{(ψ_i)* H_ij ψ_j} where 'ψ_i' is the component of the wave function on site 'i', '*' is complex conjugation, and 'H_ij' is the Hamiltonian matrix element between sites 'i' and 'j'. In Kwant this could be computed like: psi = kwant.wave_function(syst)(0)[0] # scattering state 0 from lead 0 current = np.array([2 * psi[i].conjugate() * syst.hamiltonian(i, j) * psi[j] for i, j in syst.graph]) current = current.imag In Kwant version 1.3 there will also be functionality for plotting such quantities defined on hoppings. Hope that helps, Joe
Dear Enrique, What I see in your program is that you put a potential ''edge_potential=10'' to help you to color the sites at the edge where you want to put disorder. In your case, the width is uniform. So you can achieve your goal just by putting a constraint on the coordinates of your sites: Ex: W-1>site.pos[1]>1. The most interesting case, is when the shape of your system is arbitrary. To choose the edge , you do: for site in sys.expand(graphene.shape(rec, (0, 1))): if sys.degree(site)<9 : sys[site]=edge_potential Please, notice that I used <9 and not <=9 as you did in your program (this is why everything is green in your plot). To get rid of the sites at the edge, in front of the leads, just do, after attaching the leads : for site in sys.leads[0].interface+sys.leads[1].interface: sys[site]=0 This way, your program will work for arbitrary shape, and you can think how to put disorder after this. I hope that this helps. Adel On Wed, Feb 8, 2017 at 5:31 AM, enrique colomes <kc3404@hotmail.com> wrote:
Hello,
My name is Enrique and first of all, thanks for the developping of such a useful tool.
My question is related to that one of the edge potentials. I am trying to reproduce the Haldane model with disorder in a rectangular scattering region. My goal is to introduce edge disorder as well as compute the local current density.
I have two problems:
a) Regarding the edge disorder, I cannot remove the disorder in the edges with the leads (I saw your recent previous discussion)
b) Is it possible already to compute directly from Kwant the local current density? If yes, how?
Thanks again,
Enrique
from __future__ import division # so that 1/2 == 0.5, and not 0
from scipy.sparse import spdiags from math import pi, sqrt, tanh, e import numpy as np import kwant import random import math
# For computing eigenvalues import scipy.sparse.linalg as sla
# For plotting from matplotlib import pyplot
# Define the graphene lattice graphene = kwant.lattice.honeycomb() a, b = graphene.sublattices
edge_potential=10
def make_system(L=10, W=10, pot=0.001):
t=1 tt= 0.04 * t phi=1*pi/2 edge=1 # 1 if there is edge disorder bulk=0 # 1 if there is bulk disorder nnn_hoppinga=tt*e**(1j*phi) nnn_hoppingb=tt*e**(1j*phi)
#### Define the scattering region. #### def rec(pos): x, y = pos return 0 < x < L and 0 < y < W
sys = kwant.Builder()
# w: width and pot: potential maximum of the p-n junction def potential(site): (x, y) = site.pos d = y * 2 + x * 3 return 0
sys[graphene.shape(rec, (1, 1))] = potential
# specify the hoppings of the graphene lattice nn_hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b)) nnn_hoppings_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a)) nnn_hoppings_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b)) sys[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb
# Modify the scattering region if edge==1: # change the values of the potential at the edge for site in sys.expand(graphene.shape(rec, (0, 1))): if sys.degree(site)<=9 and abs(site.pos[0])!=L-4: #abs(site.pos[0])!=19 execludes the interfaces with the leads sys[site]=edge_potential
elif bulk==1: # change the values of the potential at the edge for site in sys.expand(graphene.shape(rec, (0, 1))): if sys.degree(site)==9 and abs(site.pos[0])!=19: #abs(site.pos[0])!=19 execludes the interfaces with the leads sys[site]=edge_potential
#### Define the leads. #### # left lead sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0)))
sym0.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)]) sym0.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)])
def lead0_shape(pos): x, y = pos return (0 < y < W)
lead0 = kwant.Builder(sym0) lead0[graphene.shape(lead0_shape, (0, 0))] = -pot lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb
# right lead sym1 = kwant.TranslationalSymmetry(graphene.vec((1, 0)))
sym1.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)]) sym1.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)]) def lead1_shape(pos): x, y = pos return (0 < y < W)
lead1 = kwant.Builder(sym1) lead1[graphene.shape(lead1_shape, (0, 0))] = pot lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb
sys.attach_lead(lead1) sys.attach_lead(lead0)
return sys
# Call the main function if the script gets executed (as opposed to imported).
if __name__ == '__main__': main()t.show()
def main():
sys = make_system()
# Check that the system looks as intended.
def family_color(site): #print(sys[site]) if sys[site]==edge_potential: return 'green' else: return 'black'
def site_size(site): if sys[site]==edge_potential: return 0.35 else:return 0.2 kwant.plot(sys,site_color=family_color,site_size=site_size) # Finalize the system. sys = sys.finalized()
# Compute number of neighbors i=sys.sites.index(graphene(5,6)) all_the_neighbors=sys.graph.out_neighbors(i)
# Call the main function if the script gets executed (as opposed to imported). # See <http://docs.python.org/library/__main__.html> <http://docs.python.org/library/__main__.html>. if __name__ == '__main__': main()
-- Abbout Adel
Hello, Thanks both of you for the quick answers. I have Phyton3. It was an old line from another code. a) Clearly I didn't explain myself correctly about the edges and leads. But everything is solved now, following Adel's instructions. b) When I asked about computing it directly from Kwant it was because I read in another thread that it was going to be in Kwant1.3, but I was not sure if it was already avalaible. I will "fight" then with the code and will try to use your guidelines Joseph. Definetively doing it with neighbors will be much faster than itering the whole system. Will this new version be ready soon? :-) I will let you know my progress and will attach the code if successful. Thanks again, Enrique El 08/02/2017 a las 5:19, Abbout Adel escribió: Dear Enrique, What I see in your program is that you put a potential ''edge_potential=10'' to help you to color the sites at the edge where you want to put disorder. In your case, the width is uniform. So you can achieve your goal just by putting a constraint on the coordinates of your sites: Ex: W-1>site.pos[1]>1. The most interesting case, is when the shape of your system is arbitrary. To choose the edge , you do: for site in sys.expand(graphene.shape(rec, (0, 1))): if sys.degree(site)<9 : sys[site]=edge_potential Please, notice that I used <9 and not <=9 as you did in your program (this is why everything is green in your plot). To get rid of the sites at the edge, in front of the leads, just do, after attaching the leads : for site in sys.leads[0].interface+sys.leads[1].interface: sys[site]=0 This way, your program will work for arbitrary shape, and you can think how to put disorder after this. I hope that this helps. Adel On Wed, Feb 8, 2017 at 5:31 AM, enrique colomes <kc3404@hotmail.com<mailto:kc3404@hotmail.com>> wrote: Hello, My name is Enrique and first of all, thanks for the developping of such a useful tool. My question is related to that one of the edge potentials. I am trying to reproduce the Haldane model with disorder in a rectangular scattering region. My goal is to introduce edge disorder as well as compute the local current density. I have two problems: a) Regarding the edge disorder, I cannot remove the disorder in the edges with the leads (I saw your recent previous discussion) b) Is it possible already to compute directly from Kwant the local current density? If yes, how? Thanks again, Enrique from __future__ import division # so that 1/2 == 0.5, and not 0 from scipy.sparse import spdiags from math import pi, sqrt, tanh, e import numpy as np import kwant import random import math # For computing eigenvalues import scipy.sparse.linalg as sla # For plotting from matplotlib import pyplot # Define the graphene lattice graphene = kwant.lattice.honeycomb() a, b = graphene.sublattices edge_potential=10 def make_system(L=10, W=10, pot=0.001): t=1 tt= 0.04 * t phi=1*pi/2 edge=1 # 1 if there is edge disorder bulk=0 # 1 if there is bulk disorder nnn_hoppinga=tt*e**(1j*phi) nnn_hoppingb=tt*e**(1j*phi) #### Define the scattering region. #### def rec(pos): x, y = pos return 0 < x < L and 0 < y < W sys = kwant.Builder() # w: width and pot: potential maximum of the p-n junction def potential(site): (x, y) = site.pos d = y * 2 + x * 3 return 0 sys[graphene.shape(rec, (1, 1))] = potential # specify the hoppings of the graphene lattice nn_hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b)) nnn_hoppings_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a)) nnn_hoppings_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b)) sys[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb # Modify the scattering region if edge==1: # change the values of the potential at the edge for site in sys.expand(graphene.shape(rec, (0, 1))): if sys.degree(site)<=9 and abs(site.pos[0])!=L-4: #abs(site.pos[0])!=19 execludes the interfaces with the leads sys[site]=edge_potential elif bulk==1: # change the values of the potential at the edge for site in sys.expand(graphene.shape(rec, (0, 1))): if sys.degree(site)==9 and abs(site.pos[0])!=19: #abs(site.pos[0])!=19 execludes the interfaces with the leads sys[site]=edge_potential #### Define the leads. #### # left lead sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0))) sym0.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)]) sym0.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)]) def lead0_shape(pos): x, y = pos return (0 < y < W) lead0 = kwant.Builder(sym0) lead0[graphene.shape(lead0_shape, (0, 0))] = -pot lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb # right lead sym1 = kwant.TranslationalSymmetry(graphene.vec((1, 0))) sym1.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)]) sym1.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)]) def lead1_shape(pos): x, y = pos return (0 < y < W) lead1 = kwant.Builder(sym1) lead1[graphene.shape(lead1_shape, (0, 0))] = pot lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb sys.attach_lead(lead1) sys.attach_lead(lead0) return sys # Call the main function if the script gets executed (as opposed to imported). if __name__ == '__main__': main()t.show() def main(): sys = make_system() # Check that the system looks as intended. def family_color(site): #print(sys[site]) if sys[site]==edge_potential: return 'green' else: return 'black' def site_size(site): if sys[site]==edge_potential: return 0.35 else:return 0.2 kwant.plot(sys,site_color=family_color,site_size=site_size) # Finalize the system. sys = sys.finalized() # Compute number of neighbors i=sys.sites.index(graphene(5,6)) all_the_neighbors=sys.graph.out_neighbors(i) # Call the main function if the script gets executed (as opposed to imported). # See <http://docs.python.org/library/__main__.html><http://docs.python.org/library/__main__.html>. if __name__ == '__main__': main() -- Abbout Adel
Hi again,
b) When I asked about computing it directly from Kwant it was because I read in another thread that it was going to be in Kwant1.3, but I was not sure if it was already avalaible. I will "fight" then with the code and will try to use your guidelines Joseph. Definetively doing it with neighbors will be much faster than itering the whole system. Will this new version be ready soon? :-)
Well it is already available in the 'development' version of Kwant. If you're using 'conda' to install kwant you can easily install the development version with: conda install -c kwant kwant otherwise you can follow the instructions on the website to install the development version of Kwant "from source". There are still a couple of features that we want to implement before making the 1.3 release, but I would imagine that this release would be before the end of the month. Joe
Hello, I was successful in introducing disorder in my Haldane (second nearest neighbours) system. I am using kwant1.3, I have read the new two threads [1] and [2], and I have read [3]. However, I am still unable to compute and plot the local density current for my system in the scattering region. I would really appreciate any help. I attach the code. Once I have it I will attach it, so people can actually use it and solve future problems. Regards, Enrique [1] https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001175.html<https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001175.html><https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001175.html> [2] https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html<https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html><https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html> [3] <https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html><https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html> https://kwant-project.org/doc/dev/reference/generated/kwant.operator.Current
from scipy.sparse import spdiags from math import pi, sqrt, tanh, e
import numpy as np import kwant import random import math import tinyarray import scipy.sparse.linalg as sla
from matplotlib import pyplot
# Define the graphene lattice graphene = kwant.lattice.honeycomb() a, b = graphene.sublattices L=4 W=4 t=1 tt= 0.03 * t phi=1*pi/2 def make_system(pot=0.000): nnn_hoppinga=tt*e**(1j*phi) nnn_hoppingb=tt*e**(1j*phi) #### Define the scattering region. #### def rec(pos): x, y = pos return 0 < x < L and 0 < y < W sys = kwant.Builder() # w: width and pot: potential maximum of the p-n junction def potential(site): (x, y) = site.pos d = y * 2 + x * 3 return 0 sys[graphene.shape(rec, (1, 1))] = potential # specify the hoppings of the graphene lattice nn_hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b)) nnn_hoppings_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a)) nnn_hoppings_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b)) sys[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb #### Define the leads. #### # left lead sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0))) sym0.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)]) sym0.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)]) def lead0_shape(pos): x, y = pos return (0 < y < W) lead0 = kwant.Builder(sym0) lead0[graphene.shape(lead0_shape, (0, 0))] = -pot lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb # right lead sym1 = kwant.TranslationalSymmetry(graphene.vec((1, 0))) sym1.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)]) sym1.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)]) def lead1_shape(pos): x, y = pos return (0 < y < W) lead1 = kwant.Builder(sym1) lead1[graphene.shape(lead1_shape, (0, 0))] = pot lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb sys.attach_lead(lead1) sys.attach_lead(lead0) for site in sys.leads[0].interface+sys.leads[1].interface: sys[site]=0 return sys def jcurrent(sys): J=kwant.operator.Current(sys) wfs = kwant.wave_function(sys) wf= wfs[0] j=J(wf) kwant.plotter.map(sys, j) def main(): sys = make_system() # Finalize the system. sys = sys.finalized() # Plot local current jcurrent(sys) # Call the main function if the script gets executed (as opposed to imported). # See <http://docs.python.org/library/__main__.html><http://docs.python.org/library/__main__.html>. if __name__ == '__main__': main() El 09/02/2017 a las 2:04, Joseph Weston escribió: Hi again, b) When I asked about computing it directly from Kwant it was because I read in another thread that it was going to be in Kwant1.3, but I was not sure if it was already avalaible. I will "fight" then with the code and will try to use your guidelines Joseph. Definetively doing it with neighbors will be much faster than itering the whole system. Will this new version be ready soon? :-) Well it is already available in the 'development' version of Kwant. If you're using 'conda' to install kwant you can easily install the development version with: conda install -c kwant kwant otherwise you can follow the instructions on the website to install the development version of Kwant "from source". There are still a couple of features that we want to implement before making the 1.3 release, but I would imagine that this release would be before the end of the month. Joe
Hi,
I am using kwant1.3, I have read the new two threads [1] and [2], and I have read [3]. However, I am still unable to compute and plot the local density current for my system in the scattering region.
In general you should aim to give more information about precisely what it is that is not working. Please try and phrase your problem as "I am trying to do X; I expect the output of this code to be Y, but I receive Z". Just stating that you are unable to do something does not give anyone wishing to help enough information to be able to do so. That being said, I did notice that in the following function:
def jcurrent(sys):
J=kwant.operator.Current(sys) wfs = kwant.wave_function(sys) wf= wfs[0] j=J(wf) kwant.plotter.map(sys, j)
You attempt to plot the output of a kwant operator using 'kwant.plotter.map'. This will raise an error because 'kwant.operator.Current's return the current defined on all *hoppings* of the system, whereas 'kwant.plotter.map' only knows how to plot quantities defined on all *sites* in the system. We are currently working on adding a function 'kwant.plotter.vector_map' that will be able to plot quantities defined on hoppings (such as the current), however this is presently not available in Kwant (not even the development version). I already responded to Sverre Gulbrandsen overe here [1] where I give some examples of using current operators. Hopefully that will be enough to get you started. There will be a full tutorial on operators when Kwant 1.3 is released. Happy Kwanting, Joe [1]: https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg01150.html
Hello, Definetively I didn't help at all with my previous email. I want to check that the Haldane model in a ribbon is a gapped system in the bulk with edge states. For that purpose, I have modified my jcurrent function (attached at the end) and now I obtain correctly the current through all the hopping parameters of the system. I am trying to plot (not necessarily with Kwant) the total local density current in the x and y direction. For that I want to define for each (A/B) atom the outgoing current as the vector sum of all the bond currents in each atom (each atom has 9 relevant neighbors, except the ones in the edges). Is there any way of doing that in a simple way with the new current function? Thanks a lot for your time, Enrique def jcurrent(syst): wfs = kwant.wave_function(syst, energy=0.01) J = kwant.operator.Current(syst) wf=wfs(0)[0] j = J(wf) #for (i,j), current in zip(syst.graph, j): #print('current from site {} to {} is {}'.format(syst.sites[i].tag, syst.sites[j].tag, current)) El 07/03/2017 a las 9:25, Joseph Weston escribió: Hi, I am using kwant1.3, I have read the new two threads [1] and [2], and I have read [3]. However, I am still unable to compute and plot the local density current for my system in the scattering region. In general you should aim to give more information about precisely what it is that is not working. Please try and phrase your problem as "I am trying to do X; I expect the output of this code to be Y, but I receive Z". Just stating that you are unable to do something does not give anyone wishing to help enough information to be able to do so. That being said, I did notice that in the following function: def jcurrent(sys): J=kwant.operator.Current(sys) wfs = kwant.wave_function(sys) wf= wfs[0] j=J(wf) kwant.plotter.map(sys, j) You attempt to plot the output of a kwant operator using 'kwant.plotter.map'. This will raise an error because 'kwant.operator.Current's return the current defined on all *hoppings* of the system, whereas 'kwant.plotter.map' only knows how to plot quantities defined on all *sites* in the system. We are currently working on adding a function 'kwant.plotter.vector_map' that will be able to plot quantities defined on hoppings (such as the current), however this is presently not available in Kwant (not even the development version). I already responded to Sverre Gulbrandsen overe here [1] where I give some examples of using current operators. Hopefully that will be enough to get you started. There will be a full tutorial on operators when Kwant 1.3 is released. Happy Kwanting, Joe [1]: https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg01150.html
Hello, I was successful in introducing disorder in my Haldane (second nearest neighbours) system. I am using kwant1.3, I have read the new two threads [1] and [2], and I have read [3]. However, I am still unable to compute and plot the local density current for my system in the scattering region. I would really appreciate any help. I attach the code. Once I have it I will attach it, so people can actually use it and solve future problems. Regards, Enrique [1] https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001175.html<https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001175.html><https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001175.html> [2] https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html<https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html><https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html> [3] <https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html><https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html> https://kwant-project.org/doc/dev/reference/generated/kwant.operator.Current from scipy.sparse import spdiags from math import pi, sqrt, tanh, e import numpy as np import kwant import random import math import tinyarray import scipy.sparse.linalg as sla from matplotlib import pyplot # Define the graphene lattice graphene = kwant.lattice.honeycomb() a, b = graphene.sublattices L=4 W=4 t=1 tt= 0.03 * t phi=1*pi/2 def make_system(pot=0.000): nnn_hoppinga=tt*e**(1j*phi) nnn_hoppingb=tt*e**(1j*phi) #### Define the scattering region. #### def rec(pos): x, y = pos return 0 < x < L and 0 < y < W sys = kwant.Builder() # w: width and pot: potential maximum of the p-n junction def potential(site): (x, y) = site.pos d = y * 2 + x * 3 return 0 sys[graphene.shape(rec, (1, 1))] = potential # specify the hoppings of the graphene lattice nn_hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b)) nnn_hoppings_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a)) nnn_hoppings_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b)) sys[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb #### Define the leads. #### # left lead sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0))) sym0.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)]) sym0.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)]) def lead0_shape(pos): x, y = pos return (0 < y < W) lead0 = kwant.Builder(sym0) lead0[graphene.shape(lead0_shape, (0, 0))] = -pot lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb # right lead sym1 = kwant.TranslationalSymmetry(graphene.vec((1, 0))) sym1.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)]) sym1.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)]) def lead1_shape(pos): x, y = pos return (0 < y < W) lead1 = kwant.Builder(sym1) lead1[graphene.shape(lead1_shape, (0, 0))] = pot lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb sys.attach_lead(lead1) sys.attach_lead(lead0) for site in sys.leads[0].interface+sys.leads[1].interface: sys[site]=0 return sys def jcurrent(sys): J=kwant.operator.Current(sys) wfs = kwant.wave_function(sys) wf= wfs[0] j=J(wf) kwant.plotter.map(sys, j) def main(): sys = make_system() # Finalize the system. sys = sys.finalized() # Plot local current jcurrent(sys) # Call the main function if the script gets executed (as opposed to imported). # See <http://docs.python.org/library/__main__.html><http://docs.python.org/library/__main__.html>. if __name__ == '__main__': main() El 09/02/2017 a las 2:04, Joseph Weston escribió: Hi again, b) When I asked about computing it directly from Kwant it was because I read in another thread that it was going to be in Kwant1.3, but I was not sure if it was already avalaible. I will "fight" then with the code and will try to use your guidelines Joseph. Definetively doing it with neighbors will be much faster than itering the whole system. Will this new version be ready soon? :-) Well it is already available in the 'development' version of Kwant. If you're using 'conda' to install kwant you can easily install the development version with: conda install -c kwant kwant otherwise you can follow the instructions on the website to install the development version of Kwant "from source". There are still a couple of features that we want to implement before making the 1.3 release, but I would imagine that this release would be before the end of the month. Joe
Dear Enrique, Below is a bit of code (basically, it's not mine, it was also taken from Kwant discussions), that plots the current density. Feel free to improve it. Best wishes, Sergey from scipy.sparse import spdiags def current2d(fsys,lead=0, params=None, h=None, psi=None): x, y = numpy.array([s.pos for s in fsys.sites]).T n = len(x) xd = spdiags(x, 0, n, n) yd = spdiags(y, 0, n, n) if h is None: h = fsys.hamiltonian_submatrix(sparse=True, args=params) Jx = h * xd - xd * h Jy = h * yd - yd * h if psi is None: wf = kwant.wave_function(fsys, 0, args=params) # psi = numpy.vstack([wf(l) for l in leads]) psi = wf(lead) jx = numpy.zeros_like(psi, dtype=float) jy = numpy.zeros_like(psi, dtype=float) j2 = numpy.zeros(len(fsys.sites), dtype=float) for f in xrange(len(psi)): psi_d_conj = spdiags(psi[f].conj(), 0, n, n) jx[f] = -2*numpy.imag( psi_d_conj * (Jx * psi[f]) ) jy[f] = -2*numpy.imag( psi_d_conj * (Jy * psi[f]) ) ## plot the wave-function for the first mode in the given lead # kwant.plotter.map(fsys, numpy.absolute(psi[f]), num_lead_cells=10,a=1,vmin=0,show=False) jxtot=numpy.sum(jx,axis=0) jytot=numpy.sum(jy,axis=0) j2=numpy.sqrt(jxtot**2+jytot**2) ###add currents from all the modes and then compute the square kwant.plotter.map(fsys, j2, num_lead_cells=10,a=1,vmin=0, vmax=0.1,show=False) return jxtot , jytot # If we want to plot the current density a bit of coarse graining helps. from scipy.interpolate import griddata from scipy.ndimage.filters import gaussian_filter def smooth(x, y, jx, jy, L, W, sigma=0, downsample=1, dpu=3): x_grid = numpy.linspace(-L-1., L, 2*L*dpu+1) y_grid = numpy.linspace(-1, W+1., W*dpu+1) grid = numpy.array([_.flatten() for _ in numpy.meshgrid(x_grid,y_grid)]).T jx_grid = griddata((x,y), jx, grid, method='nearest') jx_grid = jx_grid.reshape(len(y_grid), len(x_grid)) dd = dpu*downsample #x_gr = x_grid[::dd] jx_gr = gaussian_filter(jx_grid, sigma*dpu)[::dd,::dd] jy_grid = griddata((x,y), jy, grid, method='nearest') jy_grid = jy_grid.reshape(len(y_grid), len(x_grid)) #y_gr = y_grid[::dd] jy_gr = gaussian_filter(jy_grid, sigma*dpu)[::dd,::dd] ### takes each third element return jx_gr, jy_gr # the (coarse-grained) current density looks good in a stream-plot: def plot_stream(jx, jy, ttl='J', cmap='rainbow', vmax=0.03, Wl=40): fig, ax = pyplot.subplots(figsize=(11, 7)) jval = numpy.sqrt(jx**2 + jy**2) sy, sx = jx.shape # xext = # yext = (sy-1) print (sx,sy) xgr = numpy.linspace(xmin, xmax, sx) ygr = numpy.linspace(0, Ly, sy) # jval_th = numpy.tanh(jval/vmax) ax.streamplot(xgr, ygr, jx, jy, density=2, arrowsize=1, color = numpy.tanh(jx/vmax), cmap=cmap, norm = colors.Normalize(-1,1), linewidth=numpy.tanh(jval/vmax)*2) ax.set_title(ttl) ax.margins(0.05); ax.set_aspect(1); # I use bars to indicate the position of the leads around the scattering # region, you need to adjust this to your geometry ## ax.broken_barh([(0, -5)] , (0, xext), edgecolor='grey', facecolor='grey') ## ax.broken_barh([(yext, 5)] , (0, xext), edgecolor='grey', facecolor='grey') ## ax.broken_barh([(0, Wl)] , (0, -5), edgecolor='grey', facecolor='grey') ## ax.broken_barh([(yext, -Wl)] , (0, -5), edgecolor='grey', facecolor='grey') ## ax.broken_barh([(0, yext)] , (0, xext), edgecolor='grey', facecolor='None') On 07/03/17 16:20, enrique colomes wrote:
Hello,
I was successful in introducing disorder in my Haldane (second nearest neighbours) system.
I am using kwant1.3, I have read the new two threads [1] and [2], and I have read [3]. However, I am still unable to compute and plot the local density current for my system in the scattering region.
I would really appreciate any help. I attach the code.
Once I have it I will attach it, so people can actually use it and solve future problems.
Regards,
Enrique
[3] <https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-March/001181.html> https://kwant-project.org/doc/dev/reference/generated/kwant.operator.Current
from scipy.sparse import spdiags from math import pi, sqrt, tanh, e import numpy as np import kwant import random import math import tinyarray import scipy.sparse.linalg as sla from matplotlib import pyplot # Define the graphene lattice graphene = kwant.lattice.honeycomb() a, b = graphene.sublattices
L=4 W=4 t=1 tt= 0.03 * t phi=1*pi/2
def make_system(pot=0.000):
nnn_hoppinga=tt*e**(1j*phi) nnn_hoppingb=tt*e**(1j*phi)
#### Define the scattering region. #### def rec(pos): x, y = pos return 0 < x < L and 0 < y < W
sys = kwant.Builder()
# w: width and pot: potential maximum of the p-n junction def potential(site): (x, y) = site.pos d = y * 2 + x * 3 return 0
sys[graphene.shape(rec, (1, 1))] = potential
# specify the hoppings of the graphene lattice nn_hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b)) nnn_hoppings_a = (((-1, 0), a, a), ((0, 1), a, a), ((1, -1), a, a)) nnn_hoppings_b = (((1, 0), b, b), ((0, -1), b, b), ((-1, 1), b, b)) sys[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga sys[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb
#### Define the leads. #### # left lead sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0)))
sym0.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)]) sym0.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)])
def lead0_shape(pos): x, y = pos return (0 < y < W)
lead0 = kwant.Builder(sym0) lead0[graphene.shape(lead0_shape, (0, 0))] = -pot lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga lead0[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb
# right lead sym1 = kwant.TranslationalSymmetry(graphene.vec((1, 0)))
sym1.add_site_family(graphene.sublattices[0], other_vectors=[(-1, 2)]) sym1.add_site_family(graphene.sublattices[1], other_vectors=[(-1, 2)]) def lead1_shape(pos): x, y = pos return (0 < y < W)
lead1 = kwant.Builder(sym1) lead1[graphene.shape(lead1_shape, (0, 0))] = pot lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nn_hoppings]] = -t lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_a]] = -nnn_hoppinga lead1[[kwant.builder.HoppingKind(*hopping) for hopping in nnn_hoppings_b]] = -nnn_hoppingb
sys.attach_lead(lead1) sys.attach_lead(lead0)
for site in sys.leads[0].interface+sys.leads[1].interface: sys[site]=0
return sys
def jcurrent(sys): J=kwant.operator.Current(sys) wfs = kwant.wave_function(sys) wf= wfs[0] j=J(wf) kwant.plotter.map(sys, j)
def main():
sys = make_system()
# Finalize the system. sys = sys.finalized()
# Plot local current jcurrent(sys)
# Call the main function if the script gets executed (as opposed to imported). # See<http://docs.python.org/library/__main__.html>. if __name__ == '__main__': main()
El 09/02/2017 a las 2:04, Joseph Weston escribió:
Hi again,
b) When I asked about computing it directly from Kwant it was because I read in another thread that it was going to be in Kwant1.3, but I was not sure if it was already avalaible. I will "fight" then with the code and will try to use your guidelines Joseph. Definetively doing it with neighbors will be much faster than itering the whole system. Will this new version be ready soon? :-) Well it is already available in the 'development' version of Kwant. If you're using 'conda' to install kwant you can easily install the development version with:
conda install -c kwant kwant
otherwise you can follow the instructions on the website to install the development version of Kwant "from source".
There are still a couple of features that we want to implement before making the 1.3 release, but I would imagine that this release would be before the end of the month.
Joe
participants (4)
-
Abbout Adel -
enrique colomes -
Joseph Weston -
Sergey Slizovskiy