Hi Anton, Had a go with getting the quiver plot working and from what I get it seems I've produced a graph that makes sense. My only thing I'm not 100% sure I've done right is the handling of the negative values; I assume that if I have site a and site b with a negative current value it means current is flowing from b to a, whereas if the current value is positive current is flowing from a to b. In order to account for this in the graph I plot each quiver point separately, see if the current is negative and if it is change the arrow to point into the site. I account for the magnitude of the current by setting the line width to each point to be proportional to the absolute value of the current. The attached code is below: #Plot current over hoppings wf = kwant.wave_function(syst, E,args=[gate]) print('Calculated Wavefunctions') J_0 = kwant.operator.Current(syst) psi = wf(0)[0] current = J_0(psi) orderedcurrent = [(syst.pos(i), syst.pos(j), k) for (i, j), k in zip(syst.graph, current)] #From ordered current we know that the first array contains the start location of the hopping #and the second array contains the end location thus we can get XYZ,UVW from this. X = [orderedcurrent[i][0][0] for i in range(len(orderedcurrent))] Y = [orderedcurrent[i][0][1] for i in range(len(orderedcurrent))] Z = [orderedcurrent[i][0][2] for i in range(len(orderedcurrent))] U = [orderedcurrent[i][1][0] for i in range(len(orderedcurrent))] V = [orderedcurrent[i][1][1] for i in range(len(orderedcurrent))] W = [orderedcurrent[i][1][2] for i in range(len(orderedcurrent))] #Get the current magnitude at each hopping Mag = [orderedcurrent[i][2] for i in range(len(orderedcurrent))] fig = pyplot.figure(figsize=(12,10),dpi=100) ax = fig.gca(projection='3d') #Plot each point individually, if the current magnitude is negative, plot the arrow head going into the site rather than out for i in range(len(orderedcurrent)): if Mag[i] < 0: ax.quiver(X[i], Y[i], Z[i], U[i], V[i], W[i], length = (abs(Mag[i])), pivot="tip") else: ax.quiver(X[i], Y[i], Z[i], U[i], V[i], W[i], length = (abs(Mag[i])), pivot="tail") pyplot.show() Thank You, Mitchell Greenberg ________________________________ From: anton.akhmerov@gmail.com on behalf of Anton Akhmerov Sent: Tuesday, October 3, 2017 11:22:26 PM To: Mitchell Greenberg Cc: kwant-discuss@kwant-project.org Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System Hi Mitchell, There's nothing readily available, but it should certainly be doable. Firstly you need to calculate the coordinates of all sites, combine this with the graph information, and feed all this data to e.g. the matplotlib's quiver plot ( https://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html#mpl_toolkits.mplot... ). See the recent reply of Joe ( https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg01425.html ) for assemblying the coordinates of all hoppings. Please share your code and your experience with this kind of plot if you succeed—do you think it would be universally useful? Best, Anton On Tue, Oct 3, 2017 at 3:15 PM, mitchell.greenberg mailto:mitchell.greenberg@my.jcu.edu.au> wrote: Hi Joseph and Anton, Today I had a talk with my supervisor and he wondered if it was possible to plot the raw current data before interpolation, just to see how current travels through each of the hoppings.

From previous questions I know that the current array returned from J0 (psi) is in the same from as the system graph, however I was wondering if there was a simple way to plot this array to see the current vectors at each site position.

Sorry if this question is naive. Thank You, Mitchell Greenberg Sent from Samsung Mobile -------- Original message -------- From: "mitchell.greenberg" Date:17/09/2017 7:32 PM (GMT+10:00) To: Anton Akhmerov Cc: kwant-discuss@kwant-project.orgmailto:kwant-discuss@kwant-project.org Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System This sender failed our fraud detection checks and may not be who they appear to be. Learn about spoofinghttp://aka.ms/LearnAboutSpoofing Feedbackhttp://aka.ms/SafetyTipsFeedback Hi Anton, Thanks for that. While my lattice is polyatomic I just switched site.pos to site.tag in the code I posted last message to get the tags for my lattice. Doing that I see the tags for the z axis range from 0 to 9 (makes sense as I have 10 layers), which explains why 5.3 worked when I used it but not 10. My function runs now an accurately deletes the sites without any problems. Thank You, Mitchell Greenberg Sent from Samsung Mobile -------- Original message -------- From: Anton Akhmerov Date:17/09/2017 5:11 AM (GMT+10:00) To: Mitchell Greenberg Cc: Joseph Weston ,kwant-discuss@kwant-project.orgmailto:kwant-discuss@kwant-project.org Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System Hi Mitchel, Deleting and adding sites is done by their lattice coordinates (site.tag), not their real space position (site.pos); if you want to access a site close to a certain position, you can use the method Monatomic.closest (see https://kwant-project.org/doc/1/reference/generated/kwant.lattice.Monatomic#...). Best, Anton On Sat, Sep 16, 2017 at 2:34 PM, Mitchell Greenberg mailto:mitchell.greenberg@my.jcu.edu.au> wrote:

Hi Joseph and Anton,

Thank you for the advice with getting the 2d graph up and running, its now producing results that look good. Sorry that I haven't replied in over 2 weeks, uni has gotten into full swing.

I wasn't sure if I needed to start a new topic for this (if I need to let me know) but I am now trying to remove lattice points randomly from my sample. I have created a function that can remove random lattice points (it is rather messy) however I have noticed that when I try to delete sites for layers other than at z=0 and z=5.3 I get the following error.:

KeyError: Site(kwant.lattice.Monatomic([[4.3763, 0.0, 0.0], [0.0, 3.3136, 0.0], [0.0, 0.0, 5.3]], [0.0, 0.0, 0.0], '0', 1), array([50, 50, 10]))

This error seems to pop up when I try to delete a site that doesn't exist at the specified coordinates however when I use:

Positions=[syst.sites[i].pos for i in range(len(syst.sites))]

I see that sites definitely exist for z=10 (should be 10.6 but the error output seems to round). I was therefore wondering if there is something I am missing when it comes to deleting sites and if there is an easy way to delete random sites from a 3d sample.

Thank You,

Mitchell Greenberg

________________________________ From: Joseph Weston mailto:joseph.weston08@gmail.com> Sent: Wednesday, August 30, 2017 8:33:01 PM To: Mitchell Greenberg; kwant-discuss@kwant-project.orgmailto:kwant-discuss@kwant-project.org Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System

Hi Mitchell,

...

Sorry again for the reply, I've had a look at the documentation and it seems that relwidth and abswidth control the size of the arrays. What I'm stuck on though is that my system has L = 80, W = 450 and H =52, I was expecting field to return something of the form (80,450,52,3), I must apologise again, I'm pretty new to all things python so my knowledge in the more complex areas is lacking.

The input to 'interpolate_current' is a current defined on the hoppings of your system. The output is a current defined over a realspace grid. First the current defined on the hoppings is convoluted with a bump function to smooth it out, then the resulting continuous vector field is sampled on a regular grid. The docstring explains how this is done. The interpolation grid does not in general coincide with the "grid" of sites that form your system. Typically there will be several interpolation points between neighboring sites.

The 'box' returned by 'interpolate_current' allows you to calculate the realspace position of the interpolation points.

Also the modifications to 'interpolate_current' so that it functions correctly for 3D systems just landed on the most recent development version of kwant [1].

Happy Kwanting,

Joe

[1]: https://gitlab.kwant-project.org/kwant/kwant/commit/be364e7c64c316caf86869df...

Hi Anton, The plot was relatively useful yes, we used it to see what hoppings in our system were providing the majority of transport (we were using the 10 intralayer and 5 interlayer tight binding model for Black Phosphorus). I've submitted thesis now but in my downtime I'll have a tinker with plotting it in ipyvolume, I've had a look at the docs on it and having the XYZ and UVW already it should be a case of simply calling ipyvolume. When I get it working and looking good I'll send the code in case anyone wants to do something similar down the line. Thank You, Mitchell Greenberg ________________________________ From: anton.akhmerov@gmail.com on behalf of Anton Akhmerov Sent: Monday, October 9, 2017 12:11:14 AM To: Mitchell Greenberg Cc: kwant-discuss@kwant-project.org Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System Hi Mitchell, Looks good, thanks for sharing. Are you satisfied with the usefulness of such a plot? By the way you could get a better looking result by using any plotting library that supports true 3D graphics (matplotlib does not do actual 3D rendering). Within Python world, possible options are plotly, mayavi, or ipyvolume. Cheers, Anton On Wed, Oct 4, 2017 at 10:23 AM, Mitchell Greenberg mailto:mitchell.greenberg@my.jcu.edu.au> wrote: Hi Anton, Had a go with getting the quiver plot working and from what I get it seems I've produced a graph that makes sense. My only thing I'm not 100% sure I've done right is the handling of the negative values; I assume that if I have site a and site b with a negative current value it means current is flowing from b to a, whereas if the current value is positive current is flowing from a to b. In order to account for this in the graph I plot each quiver point separately, see if the current is negative and if it is change the arrow to point into the site. I account for the magnitude of the current by setting the line width to each point to be proportional to the absolute value of the current. The attached code is below: #Plot current over hoppings wf = kwant.wave_function(syst, E,args=[gate]) print('Calculated Wavefunctions') J_0 = kwant.operator.Current(syst) psi = wf(0)[0] current = J_0(psi) orderedcurrent = [(syst.pos(i), syst.pos(j), k) for (i, j), k in zip(syst.graph, current)] #From ordered current we know that the first array contains the start location of the hopping #and the second array contains the end location thus we can get XYZ,UVW from this. X = [orderedcurrent[i][0][0] for i in range(len(orderedcurrent))] Y = [orderedcurrent[i][0][1] for i in range(len(orderedcurrent))] Z = [orderedcurrent[i][0][2] for i in range(len(orderedcurrent))] U = [orderedcurrent[i][1][0] for i in range(len(orderedcurrent))] V = [orderedcurrent[i][1][1] for i in range(len(orderedcurrent))] W = [orderedcurrent[i][1][2] for i in range(len(orderedcurrent))] #Get the current magnitude at each hopping Mag = [orderedcurrent[i][2] for i in range(len(orderedcurrent))] fig = pyplot.figure(figsize=(12,10),dpi=100) ax = fig.gca(projection='3d') #Plot each point individually, if the current magnitude is negative, plot the arrow head going into the site rather than out for i in range(len(orderedcurrent)): if Mag[i] < 0: ax.quiver(X[i], Y[i], Z[i], U[i], V[i], W[i], length = (abs(Mag[i])), pivot="tip") else: ax.quiver(X[i], Y[i], Z[i], U[i], V[i], W[i], length = (abs(Mag[i])), pivot="tail") pyplot.show() Thank You, Mitchell Greenberg ________________________________ From: anton.akhmerov@gmail.commailto:anton.akhmerov@gmail.com mailto:anton.akhmerov@gmail.com> on behalf of Anton Akhmerov mailto:anton.akhmerov%2Bkd@gmail.com> Sent: Tuesday, October 3, 2017 11:22:26 PM To: Mitchell Greenberg Cc: kwant-discuss@kwant-project.orgmailto:kwant-discuss@kwant-project.org Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System Hi Mitchell, There's nothing readily available, but it should certainly be doable. Firstly you need to calculate the coordinates of all sites, combine this with the graph information, and feed all this data to e.g. the matplotlib's quiver plot ( https://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html#mpl_toolkits.mplot... ). See the recent reply of Joe ( https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg01425.html ) for assemblying the coordinates of all hoppings. Please share your code and your experience with this kind of plot if you succeed—do you think it would be universally useful? Best, Anton On Tue, Oct 3, 2017 at 3:15 PM, mitchell.greenberg mailto:mitchell.greenberg@my.jcu.edu.au> wrote: Hi Joseph and Anton, Today I had a talk with my supervisor and he wondered if it was possible to plot the raw current data before interpolation, just to see how current travels through each of the hoppings.

From previous questions I know that the current array returned from J0 (psi) is in the same from as the system graph, however I was wondering if there was a simple way to plot this array to see the current vectors at each site position.

Sorry if this question is naive. Thank You, Mitchell Greenberg Sent from Samsung Mobile -------- Original message -------- From: "mitchell.greenberg" Date:17/09/2017 7:32 PM (GMT+10:00) To: Anton Akhmerov Cc: kwant-discuss@kwant-project.orgmailto:kwant-discuss@kwant-project.org Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System This sender failed our fraud detection checks and may not be who they appear to be. Learn about spoofinghttp://aka.ms/LearnAboutSpoofing Feedbackhttp://aka.ms/SafetyTipsFeedback Hi Anton, Thanks for that. While my lattice is polyatomic I just switched site.pos to site.tag in the code I posted last message to get the tags for my lattice. Doing that I see the tags for the z axis range from 0 to 9 (makes sense as I have 10 layers), which explains why 5.3 worked when I used it but not 10. My function runs now an accurately deletes the sites without any problems. Thank You, Mitchell Greenberg Sent from Samsung Mobile -------- Original message -------- From: Anton Akhmerov Date:17/09/2017 5:11 AM (GMT+10:00) To: Mitchell Greenberg Cc: Joseph Weston ,kwant-discuss@kwant-project.orgmailto:kwant-discuss@kwant-project.org Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System Hi Mitchel, Deleting and adding sites is done by their lattice coordinates (site.tag), not their real space position (site.pos); if you want to access a site close to a certain position, you can use the method Monatomic.closest (see https://kwant-project.org/doc/1/reference/generated/kwant.lattice.Monatomic#...). Best, Anton On Sat, Sep 16, 2017 at 2:34 PM, Mitchell Greenberg mailto:mitchell.greenberg@my.jcu.edu.au> wrote:

Hi Joseph and Anton,

Thank you for the advice with getting the 2d graph up and running, its now producing results that look good. Sorry that I haven't replied in over 2 weeks, uni has gotten into full swing.

I wasn't sure if I needed to start a new topic for this (if I need to let me know) but I am now trying to remove lattice points randomly from my sample. I have created a function that can remove random lattice points (it is rather messy) however I have noticed that when I try to delete sites for layers other than at z=0 and z=5.3 I get the following error.:

KeyError: Site(kwant.lattice.Monatomic([[4.3763, 0.0, 0.0], [0.0, 3.3136, 0.0], [0.0, 0.0, 5.3]], [0.0, 0.0, 0.0], '0', 1), array([50, 50, 10]))

This error seems to pop up when I try to delete a site that doesn't exist at the specified coordinates however when I use:

Positions=[syst.sites[i].pos for i in range(len(syst.sites))]

I see that sites definitely exist for z=10 (should be 10.6 but the error output seems to round). I was therefore wondering if there is something I am missing when it comes to deleting sites and if there is an easy way to delete random sites from a 3d sample.

Thank You,

Mitchell Greenberg

________________________________ From: Joseph Weston mailto:joseph.weston08@gmail.com> Sent: Wednesday, August 30, 2017 8:33:01 PM To: Mitchell Greenberg; kwant-discuss@kwant-project.orgmailto:kwant-discuss@kwant-project.org Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System

Hi Mitchell,

...

Sorry again for the reply, I've had a look at the documentation and it seems that relwidth and abswidth control the size of the arrays. What I'm stuck on though is that my system has L = 80, W = 450 and H =52, I was expecting field to return something of the form (80,450,52,3), I must apologise again, I'm pretty new to all things python so my knowledge in the more complex areas is lacking.

The input to 'interpolate_current' is a current defined on the hoppings of your system. The output is a current defined over a realspace grid. First the current defined on the hoppings is convoluted with a bump function to smooth it out, then the resulting continuous vector field is sampled on a regular grid. The docstring explains how this is done. The interpolation grid does not in general coincide with the "grid" of sites that form your system. Typically there will be several interpolation points between neighboring sites.

The 'box' returned by 'interpolate_current' allows you to calculate the realspace position of the interpolation points.

Also the modifications to 'interpolate_current' so that it functions correctly for 3D systems just landed on the most recent development version of kwant [1].

Happy Kwanting,

Joe

[1]: https://gitlab.kwant-project.org/kwant/kwant/commit/be364e7c64c316caf86869df...