[Kwant] Visualize current through a 2D cut for a 3D system
Hello everyone, Is there currently way to visualize the current through a 2D cut for a 3D system that does not involve messing with the kwant code? I have a 3D system and would like to see the differences in the current when I switch some specific hoppings on and off. Regards, -- Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elch...@auth.gr - tel:+30 2310 998109
Re: [Kwant] Visualize current through a 2D cut for a 3D system
Dear Eleni, I think this should be doable with not too much of low level work. * Firstly you probably will want to select only the hoppings that are close to the 2D cut. You can do this by using "where" parameter to Current operator (see https://kwant-project.org/doc/1/reference/generated/kwant.operator.Current#kwant.operator.Current) * After computing the current across those hoppings, you would need to call kwant.plotter.interpolate_current. This function only does the interpolation, but no plotting. Having that interpolation in 3D you need to slice it and probably take a projection onto 2D plane. * Finally, the resulting 2D array you can feed to kwant.plotter.streamplot kwant.plotter.current is essentially a straightforward combination of kwant.plotter.interpolate_current and kwant.plotter.streamplot Also: your problem sounds like a useful thing to do. Please let me know if you succeed, and what you ended up doing (or if you face any further problems as well). Best, Anton On Fri, Mar 23, 2018 at 5:22 PM, wrote: > Hello everyone, > > Is there currently way to visualize the current through a 2D cut for a 3D > system that does not involve messing with the kwant code? > > I have a 3D system and would like to see the differences in the current when > I switch some specific hoppings on and off. > > > Regards, > > > > -- > Dr. Eleni Chatzikyriakou > Computational Physics lab > Aristotle University of Thessaloniki > elch...@auth.gr - tel:+30 2310 998109 >
Re: [Kwant] Visualize current through a 2D cut for a 3D system
Hello Anton, Thank you very much for your reply. I am having a hard time with the "where" argument in kwant.operator.Current. My system is not finalized. How can I find the sites associated with the hopping? By their tag? Their position? Do you have a piece of code for that? Note that I am using TBModels and Wannier90 and the lattice has one sublattice for each WF (or site, or orbital). norb is not defined, but I see that for x in lattice.sublattices: x.norbs=1 works and looks correct. Regards, Eleni Quoting Anton Akhmerov : Dear Eleni, I think this should be doable with not too much of low level work. * Firstly you probably will want to select only the hoppings that are close to the 2D cut. You can do this by using "where" parameter to Current operator (see https://kwant-project.org/doc/1/reference/generated/kwant.operator.Current#kwant.operator.Current) * After computing the current across those hoppings, you would need to call kwant.plotter.interpolate_current. This function only does the interpolation, but no plotting. Having that interpolation in 3D you need to slice it and probably take a projection onto 2D plane. * Finally, the resulting 2D array you can feed to kwant.plotter.streamplot kwant.plotter.current is essentially a straightforward combination of kwant.plotter.interpolate_current and kwant.plotter.streamplot Also: your problem sounds like a useful thing to do. Please let me know if you succeed, and what you ended up doing (or if you face any further problems as well). Best, Anton On Fri, Mar 23, 2018 at 5:22 PM, wrote: Hello everyone, Is there currently way to visualize the current through a 2D cut for a 3D system that does not involve messing with the kwant code? I have a 3D system and would like to see the differences in the current when I switch some specific hoppings on and off. Regards, -- Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elch...@auth.gr - tel:+30 2310 998109 -- Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elch...@auth.gr - tel:+30 2310 998109
Re: [Kwant] Visualize current through a 2D cut for a 3D system
Hello everyone, I'm not sure what I am doing wrong, but here is the output when using interpolate_current(): kwant_sys = wraparound.wraparound(kwant_model).finalized() ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False) evecs = sla.eigsh(ham_mat, k=48, which='SM')[1] J = kwant.operator.Current(kwant_sys) current = J(evecs[:, 9]) IJ = kwant.plotter.interpolate_current(kwant_sys, J) Traceback (most recent call last): File "", line 1, in File "/apps/applications/python/anaconda3/5.0.1/lib/python3.6/site-packages/kwant/plotter.py", line 1852, in interpolate_current if len(current) != syst.graph.num_edges: TypeError: object of type 'kwant.operator.Current' has no len() current array([ 1.04495760e-10, 1.04407946e-10, 4.08360326e-10, ..., -5.89810275e-10, -2.77384947e-10, 1.41194501e-10]) I was also looking at this thread: https://mailman-mail5.webfaction.com/pipermail/kwant-discuss/2017-October/001466.html and I think I will try also using ipyvolume or mayavi. But I would like to know if there is a problem with this being a closed system. Regards, Eleni Quoting Anton Akhmerov : Dear Eleni, I think this should be doable with not too much of low level work. * Firstly you probably will want to select only the hoppings that are close to the 2D cut. You can do this by using "where" parameter to Current operator (see https://kwant-project.org/doc/1/reference/generated/kwant.operator.Current#kwant.operator.Current) * After computing the current across those hoppings, you would need to call kwant.plotter.interpolate_current. This function only does the interpolation, but no plotting. Having that interpolation in 3D you need to slice it and probably take a projection onto 2D plane. * Finally, the resulting 2D array you can feed to kwant.plotter.streamplot kwant.plotter.current is essentially a straightforward combination of kwant.plotter.interpolate_current and kwant.plotter.streamplot Also: your problem sounds like a useful thing to do. Please let me know if you succeed, and what you ended up doing (or if you face any further problems as well). Best, Anton On Fri, Mar 23, 2018 at 5:22 PM, wrote: Hello everyone, Is there currently way to visualize the current through a 2D cut for a 3D system that does not involve messing with the kwant code? I have a 3D system and would like to see the differences in the current when I switch some specific hoppings on and off. Regards, -- Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elch...@auth.gr - tel:+30 2310 998109 -- Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elch...@auth.gr - tel:+30 2310 998109
Re: [Kwant] Visualize current through a 2D cut for a 3D system
Hi, On 04/10/2018 01:50 PM, elch...@auth.gr wrote: > Hello everyone, > > I'm not sure what I am doing wrong, but here is the output when using > interpolate_current(): > > kwant_sys = wraparound.wraparound(kwant_model).finalized() ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False) evecs = sla.eigsh(ham_mat, k=48, which='SM')[1] J = kwant.operator.Current(kwant_sys) current = J(evecs[:, 9]) IJ = kwant.plotter.interpolate_current(kwant_sys, J) You're passing "J" to "interpolate current", when you should pass "current" > > and I think I will try also using ipyvolume or mayavi. But I would > like to know if there is a problem with this being a closed system. No, no problem. Of course if you don't break time reversal symmetry (e.g. with magnetic field) then you're eigenstates will have 0 current everywhere. Happy Kwanting, Joe
Re: [Kwant] Visualize current through a 2D cut for a 3D system
Hello Joseph,
Sorry I copied the wrong pieces of code:
ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False)
ev = sla.eigsh(ham_mat, k=48, which='SM')
evecs = ev[1]
lat1 = lattice.sublattices[1]
lat3 = lattice.sublattices[3]
J = kwant.operator.Current(kwant_sys, where=[(lat1(1, 0, 0), lat3(0, 0, 0))])
current = J(evecs[:, 1])
IJ = kwant.plotter.interpolate_current(kwant_sys, current)
---
ValueErrorTraceback (most recent call last)
in ()
5 J = kwant.operator.Current(kwant_sys, where=[(lat1(1, 0, 0),
lat3(0, 0, 0))])
6 current = J(evecs[:, 1])
> 7 IJ = kwant.plotter.interpolate_current(kwant_sys, current)
8 print(current)
9
~\Anaconda3\lib\site-packages\kwant\plotter.py in
interpolate_current(syst, current, relwidth, abswidth, n)
1851
1852 if len(current) != syst.graph.num_edges:
-> 1853 raise ValueError("Current and hoppings arrays do not
have the same"
1854 " length.")
1855
ValueError: Current and hoppings arrays do not have the same length.
print(current)
[7.65418274e-10]
--
---
ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False)
ev = sla.eigsh(ham_mat, k=48, which='SM')
evecs = ev[1]
J = kwant.operator.Current(kwant_sys)
current = J(evecs[:, 1])
IJ = kwant.plotter.interpolate_current(kwant_sys, current)
---
ValueErrorTraceback (most recent call last)
in ()
6 J = kwant.operator.Current(kwant_sys)
7 current = J(evecs[:, 1])
> 8 IJ = kwant.plotter.interpolate_current(kwant_sys, current)
9
~\Anaconda3\lib\site-packages\kwant\plotter.py in
interpolate_current(syst, current, relwidth, abswidth, n)
1900 # this check. This check is done here to keep changes local
1901 if dim != 2:
-> 1902 raise ValueError("'interpolate_current' only works for
2D systems.")
1903 factor = (3 / np.pi) / (width / 2)
1904 scale = 2 / width
ValueError: 'interpolate_current' only works for 2D systems.
---
I was able to do a points3d() plot with the 3D data in Mayavi using
the following piece of code:
---
ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False)
ev = sla.eigsh(ham_mat, k=48, which='SM')
evecs = ev[1]
prob_dens = np.abs(evecs[:, 47])**2
x = np.array([])
y = np.array([])
z = np.array([])
for i in range(len(kwant_sys.sites)):
x = np.append(x,[kwant_sys.pos(i)[0]], axis=0)
y = np.append(y,[kwant_sys.pos(i)[1]], axis=0)
z = np.append(z,[kwant_sys.pos(i)[2]], axis=0)
pts = mlab.points3d(x, y, z, prob_dens)
mayavi.mlab.scalarbar(object=pts)
mayavi.mlab.xlabel("X", object = pts)
mlab.savefig("pts_47.png")
mlab.show()
However, a contour3d would be more appropriate, but it looks like that
requires regularily spaced points.
On the other hand, an unstructured grid looks a bit more cumbersome to create.
At this point, the easiest thing looks like copying from kwant code
directly into my notebook for interpolation etc.
Eleni
Quoting Joseph Weston :
Hi,
On 04/10/2018 01:50 PM, elch...@auth.gr wrote:
Hello everyone,
I'm not sure what I am doing wrong, but here is the output when using
interpolate_current():
kwant_sys = wraparound.wraparound(kwant_model).finalized()
ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False)
evecs = sla.eigsh(ham_mat, k=48, which='SM')[1]
J = kwant.operator.Current(kwant_sys)
current = J(evecs[:, 9])
IJ = kwant.plotter.interpolate_current(kwant_sys, J)
You're passing "J" to "interpolate current", when you should pass "current"
and I think I will try also using ipyvolume or mayavi. But I would
like to know if there is a problem with this being a closed system.
No, no problem. Of course if you don't break time reversal symmetry
(e.g. with magnetic field) then you're eigenstates will have 0 current
everywhere.
Happy Kwanting,
Joe
--
Dr. Eleni Chatzikyriakou
Computational Physics lab
Aristotle University of Thessaloniki
elch...@auth.gr - tel:+30 2310 998109
Re: [Kwant] Visualize current through a 2D cut for a 3D system
Hi,
> Hello Joseph,
>
> Sorry I copied the wrong pieces of code:
>
>
> ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False)
> ev = sla.eigsh(ham_mat, k=48, which='SM')
> evecs = ev[1]
> lat1 = lattice.sublattices[1]
> lat3 = lattice.sublattices[3]
> J = kwant.operator.Current(kwant_sys, where=[(lat1(1, 0, 0), lat3(0,
> 0, 0))])
> current = J(evecs[:, 1])
> IJ = kwant.plotter.interpolate_current(kwant_sys, current)
> ---
>
> ValueError Traceback (most recent call
> last)
> in ()
> 5 J = kwant.operator.Current(kwant_sys, where=[(lat1(1, 0, 0),
> lat3(0, 0, 0))])
> 6 current = J(evecs[:, 1])
> > 7 IJ = kwant.plotter.interpolate_current(kwant_sys, current)
> 8 print(current)
> 9
>
> ~\Anaconda3\lib\site-packages\kwant\plotter.py in
> interpolate_current(syst, current, relwidth, abswidth, n)
> 1851
> 1852 if len(current) != syst.graph.num_edges:
> -> 1853 raise ValueError("Current and hoppings arrays do not
> have the same"
> 1854 " length.")
> 1855
>
> ValueError: Current and hoppings arrays do not have the same length.
>
> print(current)
> [7.65418274e-10]
> --
Ah right, you are only calculating the current over a single hopping,
but 'interpolate_current' requires the current over *all* hoppings. What
extra information do you hope to get by plotting in realspace the
current across a single hopping?
This reason this restriction is in place is that 'interpolate_current'
accepts just an array of numbers for the current. This means that there
is no way for 'interpolate_current' to know which number corresponds to
the current across which hopping without further information. This
'further information' is the fact that the array of currents is in the
same order as 'syst.graph'; if we don't specify all the hoppings there's
know way to know which ones we've specified.
> ---
> ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False)
> ev = sla.eigsh(ham_mat, k=48, which='SM')
> evecs = ev[1]
> J = kwant.operator.Current(kwant_sys)
> current = J(evecs[:, 1])
> IJ = kwant.plotter.interpolate_current(kwant_sys, current)
> ---
>
> ValueError Traceback (most recent call
> last)
> in ()
> 6 J = kwant.operator.Current(kwant_sys)
> 7 current = J(evecs[:, 1])
> > 8 IJ = kwant.plotter.interpolate_current(kwant_sys, current)
> 9
>
> ~\Anaconda3\lib\site-packages\kwant\plotter.py in
> interpolate_current(syst, current, relwidth, abswidth, n)
> 1900 # this check. This check is done here to keep
> changes local
> 1901 if dim != 2:
> -> 1902 raise ValueError("'interpolate_current' only works for
> 2D systems.")
> 1903 factor = (3 / np.pi) / (width / 2)
> 1904 scale = 2 / width
>
> ValueError: 'interpolate_current' only works for 2D systems.
>
> ---
The restriction of 'interpolate_current' to 2D was lifted a while ago
(https://gitlab.kwant-project.org/kwant/kwant/merge_requests/166) but we
have not made an official release since. If you're feeling adventurous
you can install the Kwant directly from the 'master' branch of the git
repository (post back here if you don't know how).
Happy Kwanting,
Joe
Re: [Kwant] Visualize current through a 2D cut for a 3D system
Hello all,
Is it also possible to help me with this function to match first
nearest neighbor sites for the where argument in kwant.operator.Current:
def first_neighbors(site_to, site_from):
for i in site_from.family.neighbors(1):
if(??):
return True
else
return False
[...]
J = kwant.operator.Current(kwant_sys, where=first_neighbors, sum=True]
Thank you
Quoting Joseph Weston :
Hi,
Hello Joseph,
Sorry I copied the wrong pieces of code:
ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False)
ev = sla.eigsh(ham_mat, k=48, which='SM')
evecs = ev[1]
lat1 = lattice.sublattices[1]
lat3 = lattice.sublattices[3]
J = kwant.operator.Current(kwant_sys, where=[(lat1(1, 0, 0), lat3(0,
0, 0))])
current = J(evecs[:, 1])
IJ = kwant.plotter.interpolate_current(kwant_sys, current)
---
ValueError Traceback (most recent call
last)
in ()
5 J = kwant.operator.Current(kwant_sys, where=[(lat1(1, 0, 0),
lat3(0, 0, 0))])
6 current = J(evecs[:, 1])
> 7 IJ = kwant.plotter.interpolate_current(kwant_sys, current)
8 print(current)
9
~\Anaconda3\lib\site-packages\kwant\plotter.py in
interpolate_current(syst, current, relwidth, abswidth, n)
1851
1852 if len(current) != syst.graph.num_edges:
-> 1853 raise ValueError("Current and hoppings arrays do not
have the same"
1854 " length.")
1855
ValueError: Current and hoppings arrays do not have the same length.
print(current)
[7.65418274e-10]
--
Ah right, you are only calculating the current over a single hopping,
but 'interpolate_current' requires the current over *all* hoppings. What
extra information do you hope to get by plotting in realspace the
current across a single hopping?
This reason this restriction is in place is that 'interpolate_current'
accepts just an array of numbers for the current. This means that there
is no way for 'interpolate_current' to know which number corresponds to
the current across which hopping without further information. This
'further information' is the fact that the array of currents is in the
same order as 'syst.graph'; if we don't specify all the hoppings there's
know way to know which ones we've specified.
---
ham_mat = kwant_sys.hamiltonian_submatrix(sparse=False)
ev = sla.eigsh(ham_mat, k=48, which='SM')
evecs = ev[1]
J = kwant.operator.Current(kwant_sys)
current = J(evecs[:, 1])
IJ = kwant.plotter.interpolate_current(kwant_sys, current)
---
ValueError Traceback (most recent call
last)
in ()
6 J = kwant.operator.Current(kwant_sys)
7 current = J(evecs[:, 1])
> 8 IJ = kwant.plotter.interpolate_current(kwant_sys, current)
9
~\Anaconda3\lib\site-packages\kwant\plotter.py in
interpolate_current(syst, current, relwidth, abswidth, n)
1900 # this check. This check is done here to keep
changes local
1901 if dim != 2:
-> 1902 raise ValueError("'interpolate_current' only works for
2D systems.")
1903 factor = (3 / np.pi) / (width / 2)
1904 scale = 2 / width
ValueError: 'interpolate_current' only works for 2D systems.
---
The restriction of 'interpolate_current' to 2D was lifted a while ago
(https://gitlab.kwant-project.org/kwant/kwant/merge_requests/166) but we
have not made an official release since. If you're feeling adventurous
you can install the Kwant directly from the 'master' branch of the git
repository (post back here if you don't know how).
Happy Kwanting,
Joe
--
Dr. Eleni Chatzikyriakou
Computational Physics lab
Aristotle University of Thessaloniki
elch...@auth.gr - tel:+30 2310 998109
Re: [Kwant] Visualize current through a 2D cut for a 3D system
Hi, > > Is it also possible to help me with this function to match first > nearest neighbor sites for the where argument in kwant.operator.Current: > > def first_neighbors(site_to, site_from): > for i in site_from.family.neighbors(1): > if(??): > return True > else > return False > > [...] > > J = kwant.operator.Current(kwant_sys, where=first_neighbors, sum=True] Probably the following: def first_neighbors(to_site, from_site): return to_site in from_site.neighbors(1) 'neighbors' returns (per the documentation [1]) a list of sites and we can just use the 'in' operator to check for the requested site. Happy Kwanting, Joe [1]: https://kwant-project.org/doc/1.0/reference/generated/kwant.lattice.Polyatomic#kwant.lattice.Polyatomic.neighbors
Re: [Kwant] Visualize current through a 2D cut for a 3D system
Hello Joseph, One problem that I see is that I can only get neighbors() from Builders and lattices (using site.family). Using something like this: def first_neighbors(to_site, from_site): return to_site in from_site.family.neighbors(1) is a good idea, but am I really comparing two sites? for i in x.family.neighbors(1): print(i) HoppingKind((1, 0, 0), ) HoppingKind((0, 1, 0), ) --- Also, since this is very expensive computationally, is there any other way to find the hopping distance between two sites? Eleni Quoting Joseph Weston : Hi, Is it also possible to help me with this function to match first nearest neighbor sites for the where argument in kwant.operator.Current: def first_neighbors(site_to, site_from): for i in site_from.family.neighbors(1): if(??): return True else return False [...] J = kwant.operator.Current(kwant_sys, where=first_neighbors, sum=True] Probably the following: def first_neighbors(to_site, from_site): return to_site in from_site.neighbors(1) 'neighbors' returns (per the documentation [1]) a list of sites and we can just use the 'in' operator to check for the requested site. Happy Kwanting, Joe [1]: https://kwant-project.org/doc/1.0/reference/generated/kwant.lattice.Polyatomic#kwant.lattice.Polyatomic.neighbors -- Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elch...@auth.gr - tel:+30 2310 998109
Re: [Kwant] Visualize current through a 2D cut for a 3D system
Hi again, > Hello Joseph, > > One problem that I see is that I can only get neighbors() from > Builders and lattices (using site.family). > > Using something like this: > > def first_neighbors(to_site, from_site): > return to_site in from_site.family.neighbors(1) > > is a good idea, but am I really comparing two sites? > > > for i in x.family.neighbors(1): > print(i) > > HoppingKind((1, 0, 0), ) > HoppingKind((0, 1, 0), ) > --- > Ah, sorry I was confused. > Also, since this is very expensive computationally, is there any other > way to find the hopping distance between two sites? Sure! If the sites are from the same lattice then you can just take the difference of the "tags" of the sites: delta = to_site.tag - from_site.tag 'delta' will be an integer array: the position of the site in the basis of lattice vectors. If the sites are from different lattices then you probably can't do better than comparing the realspace positions: realspace_delta = to_site.pos - from_site.pos Happy Kwanting, Joe
Re: [Kwant] Visualize current through a 2D cut for a 3D system
Thank you. That's great. Quoting Joseph Weston : Hi again, Hello Joseph, One problem that I see is that I can only get neighbors() from Builders and lattices (using site.family). Using something like this: def first_neighbors(to_site, from_site): return to_site in from_site.family.neighbors(1) is a good idea, but am I really comparing two sites? for i in x.family.neighbors(1): print(i) HoppingKind((1, 0, 0), ) HoppingKind((0, 1, 0), ) --- Ah, sorry I was confused. Also, since this is very expensive computationally, is there any other way to find the hopping distance between two sites? Sure! If the sites are from the same lattice then you can just take the difference of the "tags" of the sites: delta = to_site.tag - from_site.tag 'delta' will be an integer array: the position of the site in the basis of lattice vectors. If the sites are from different lattices then you probably can't do better than comparing the realspace positions: realspace_delta = to_site.pos - from_site.pos Happy Kwanting, Joe -- Dr. Eleni Chatzikyriakou Computational Physics lab Aristotle University of Thessaloniki elch...@auth.gr - tel:+30 2310 998109
